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^ ! Abstract 

We describe the new version (v2.38j) of the code hfodd which solves the nuclear Skyrme- 
Hartree-Fock or Skyrme-Hartree-Fock-Bogolyubov problem by using the Cartesian deformed 
harmonic-oscillator basis. In the new version, we have implemented: (i) projection on good 
angular momentum (for the Hartree-Fock states), (ii) calculation of the GCM kernels, (iii) 
calculation of matrix elements of the Yukawa interaction, (iv) the BCS solutions for state- 
dependent pairing gaps, (v) the HFB solutions for broken simplex symmetry, (vi) calculation 
of Bohr deformation parameters, (vii) constraints on the Schiff moments and scalar multipole 
moments, (viii) the Dj h transformations and rotations of wave functions, (ix) quasiparticle 
blocking for the HFB solutions in odd and odd-odd nuclei, (x) the Broyden method to accelerate 
the convergence, (xi) the Lipkin-Nogami method to treat pairing correlations, (xii) the exact 
Coulomb exchange term, (xiii) several utility options, and we have corrected two insignificant 
errors. 

PACS numbers: 07.05.T, 21.60.-n, 21.60.Jz 

NEW VERSION PROGRAM SUMMARY 

Title of the program: hfodd (v2.38j) 



1 E-mail: jacek.dobaczewski@fuw.edu.pl 



Catalogue number: 



Program obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland 
(see application form in this issue) 

Reference in CPC for earlier version of program: J. Dobaczewski and P. Olbratowski, Comput. 
Phys. Commun. 167 (2005) 214 (v2.08k). 

Catalogue number of previous version: ADFL_v2_l 

Licensing provisions: none 

Does the new version supersede the previous one: yes 

Computers on which the program has been tested: Pentium-Ill, AMD-Athlon, AMD-Opteron 

Operating systems: UNIX, LINUX, Windows xp 

Programming language used: FORTRAN- 77 and FORTRAN-90 

Memory required to execute with typical data: 10 Mwords 

No. of bits in a word: The code is written in single-precision for the use on a 64-bit processor. 
The compiler option -r8 or +autodblpad (or equivalent) has to be used to promote all real 
and complex single-precision floating-point items to double precision when the code is used on 
a 32-bit machine. 

Has the code been vectorised?: Yes 

No. of lines in distributed program: 69 085 (of which 29 602 are comments and separators) 

Keywords: Hartree-Fock; Hartree-Fock-Bogolyubov; Skyrme interaction; Self-consistent mean- 
field; Nuclear many-body problem; Superdeformation; Quadrupole deformation; Octupole de- 
formation; Pairing; Nuclear radii; Single-particle spectra; Nuclear rotation; High-spin states; 
Moments of inertia; Level crossings; Harmonic oscillator; Coulomb field; Pairing; Point sym- 
metries; Yukawa interaction; Angular-momentum projection; Generator Coordinate Method; 
Schiff moments 

Nature of physical problem 

The nuclear mean-field and an analysis of its symmetries in realistic cases are the main in- 
gredients of a description of nuclear states. Within the Local Density Approximation, or for 
a zero-range velocity-dependent Skyrme interaction, the nuclear mean-field is local and ve- 
locity dependent. The locality allows for an effective and fast solution of the self-consistent 
Hartree-Fock equations, even for heavy nuclei, and for various nucleonic (n-particle n-hole) 
configurations, deformations, excitation energies, or angular momenta. Similar Local Density 
Approximation in the particle-particle channel, which is equivalent to using a zero-range interac- 
tion, allows for a simple implementation of pairing effects within the Hartree-Fock-Bogolyubov 
method. 
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Method of solution 

The program uses the Cartesian harmonic oscillator basis to expand single-particle or single- 
quasiparticle wave functions of neutrons and protons interacting by means of the Skyrme effec- 
tive interaction and zero-range pairing interaction. The expansion coefficients are determined 
by the iterative diagonalization of the mean field Hamiltonians or Routhians which depend 
non-linearly on the local neutron and proton densities. Suitable constraints are used to obtain 
states corresponding to a given configuration, deformation or angular momentum. The method 
of solution has been presented in: J. Dobaczewski and J. Dudek, Comput. Phys. Commun. 102 
(1997) 166. 

Summary of revisions 

1. Projection on good angular momentum (for the Hartree-Fock states) has been imple- 
mented. 

2. Calculation of the GCM kernels has been implemented. 

3. Calculation of matrix elements of the Yukawa interaction has been implemented. 

4. The BCS solutions for state-dependent pairing gaps have been implemented. 

5. The HFB solutions for broken simplex symmetry have been implemented. 

6. Calculation of Bohr deformation parameters has been implemented. 

7. Constraints on the Schiff moments and scalar multipole moments have been implemented. 

8. The Djh transformations and rotations of wave functions have been implemented. 

9. The quasiparticle blocking for the HFB solutions in odd and odd-odd nuclei has been 
implemented. 

10. The Broyden method to accelerate the convergence has been implemented. 

11. The Lipkin-Nogami method to treat pairing correlations has been implemented. 

12. The exact Coulomb exchange term has been implemented. 

13. Several utility options have been implemented. 

14. Two insignificant errors have been corrected. 

Restrictions on the complexity of the problem 

The main restriction is the CPU time required for calculations of heavy deformed nuclei and 
for a given precision required. 

Typical running time 

One Hartree-Fock iteration for the superdeformed, rotating, parity conserving state of 1 ||Dyg6 
takes about six seconds on the AMD-Athlon 1600+ processor. Starting from the Woods- 
Saxon wave functions, about fifty iterations are required to obtain the energy converged within 
the precision of about 0.1 keV. In case when every value of the angular velocity is converged 
separately, the complete superdeformed band with precisely determined dynamical moments 
can be obtained within forty minutes of CPU on the AMD-Athlon 1600+ processor. This 
time can be often reduced by a factor of three when a self-consistent solution for a given 
rotational frequency is used as a starting point for a neighboring rotational frequency. 

Unusual features of the program 

The user must have access to (i) an implementation of the BLAS (Basic Linear Algebra Sub- 
routines), (ii) the NAGLIB subroutine f02axe, or LAPACK subroutines zhpev, zhpevx, or 
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zheevr, which diagonalize complex hermitian matrices, and (iii) the LINPACK subroutines 
zgedi and zgeco, which invert arbitrary complex matrices and calculate determinants, or 
provide another set of subroutines that can perform such a tasks. The LAPACK and LIN- 
PACK subroutines and an unoptimized version of the BLAS can be obtained from the Netlib 
Repository at the University of Tennessee, Knoxville: http://www.netlib.org/. 

LONG WRITE-UP 



1 Introduction 

The method of solving the Hartree-Fock (HF) equations in the Cartesian harmonic oscillator 
(HO) basis was described in the previous publication, Ref. P (I). Four versions of the code 
hfodd were previously published: (vl.60r) [2j (II), (vl.75r) [3] (III), (v2.08i) @] (IV), and 
(v2.08k) [5] (V). The User's Guide for version (v2.08k) is available in Ref. [6] and the code 
home page is at http://www.fuw.edu.pl/~dobaczew/hfodd/hfodd.html. The present paper 
is a long write-up of the new version (v2.38j) of the code hfodd. This extended version 
features the angular-momentum projection, calculations of the generator-coordinate-method 
(GCM) kernels, and several other major modifications, and is fully compatible with all previous 
versions. 

Information provided in previous publications [1] [5] remains valid, unless explicitly men- 
tioned in the present long write-up. Below we refer by Roman numbers (I)-(V) to section and 
equation numbers in these previous publications 

In Section [2] we review modifications introduced in version (v2.38j) of the code hfodd. 
Section [3] lists all additional new input keywords and data values, introduced in version (v2.38j). 
The structure of the input data file remains the same as in the previous versions, see Section II-3. 
Similarly, all previously introduced keywords and data values retain their validity and meaning. 

2 Modifications introduced in version (v2.38j) 

2.1 Projection on good angular momentum and calculation of the GCM kernels 

The code hfodd (v2.38j) contains a new option allowing for the angular-momentum projection 
(AMP) after variation of an arbitrary symmetry- unrestricted Slater determinant \§) provided 
by the code. This includes projection of cranked time-reversal symmetry breaking HF states 
which is of particular interest in high-spin applications [H [91 [TOj EH IT2"] . 

The angular-momentum conserving wave function is obtained from the Slater determinant 
|$) by applying the standard S0(3) AMP operator [T3"l [T4"] : 

\IMK) = Pl IK \<S>) = J Di* K (Q) R(n)\*) dn, (1) 

projecting onto angular momentum I with projection M and K along the laboratory and 
intrinsic z-axes, respectively. The symbol Q labels here a set of three Euler angles (a,(3,j), 
D MK (Jl) is the Wigner function and R(Q) = e~ iaIz e~ l/3Iy e~ l/yIz is the active (body) rotation 
operator [T5] . 
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The intrinsic quantum number K is, in general, not conserved. The i^-mixing is taken into 
account by assuming the following form for the eigenstates 



\i ] IM) = Y.9?k\IMK)=Y.9?kP I mk\3>)- (2) 

K K 

The mixing coefficients gf^ are determined by a minimization of energy, which amounts to 
solving the following Hill- Wheeler (H-W) eigenvalue problem: 



Hkk'9ik' = E i X! N KK'gm', (3) 

K' K 



separately for each angular momentum /. The Hamiltonian and norm matrix elements entering 
the H-W equation (j3J) equal: 

h kk , = ($\HPk K ,m = / dil D$ K ,{Sl) W(O), (4) 
(®\Pkk>\®) = / dSl D% K ,(n) AT(Q), (5) 



N. 



KK' ~ 

where 

H{$1) = ($\HR(n)\$), (6) 
AT(n) = ($1^)1$), (7) 

denote the Hamiltonian and norm kernels, respectively. 

The H-W eigenvalue problem ([3]) should be handled with care due to its overcompleteness. 
This difficulty is overcome in the code by solving the problem ([3]) in the collective basis spanned 
by the natural states: 

\IM)W = — *— YvPlIMK). (8) 
These states are constructed of the eigenstates of the norm matrix: 

n kk'Vip = n m ffp , (9) 

K> 

having non-zero n m > eigenvalues. More precisely, due to numerical stability reasons, the 
collective subspace is constructed by using only n = l,2,...,m max eigenstates having n m > 
(, where ( is an externally provided basis cut-off parameter. Final diagonalization of the 
Hamiltonian matrix is performed in the m max -dimensional collective subspace defined as: 

"•max 

\i-IM)= £ /«|/M) (m) . (10) 

m=l 

The code provides, for each angular momentum, eigenenergies and mixing coefficients recalcu- 
lated to the representation defined by Eq. (j2J) according to the following formula: 

mmax f (i) rl ( m ) 

« = £ (11) 
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The cornerstone of the AMP scheme described above is a calculation of the Hamiltonian 
([!)]) and norm (0) kernels. A prerequisite for this calculation is a spatial rotation of the Slater 
determinant. It amounts to rotating independently each single-particle state R(Q) ipi(f,a) = 
<Pi(f,a): 

y5i(i) 

£i(2) 



R(n)\$) = |$) 



1 



.4! 



02(2) 



0a(1) 
<Pa(2) 



0a(A) 
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This particular task is performed by the subroutine rotwav. Taking advantage of a fact 
that the hfodd is coded in the Cartesian HO basis^l ipn x (x)ipn y (y)?Pn z (z) , cf. Eq. (1-71), the 
three-dimensional rotation is factorized into three independent and numerically equivalent one- 
dimensional rotations around z—, y— and again z— axis, respectively. To be more specific, 
in the first step, each single-particle wave function is rotated around the z-axis by the Euler 
angle 7, and the resulting rotated wave function is expanded in the original Cartesian harmonic 
oscillator basis: 



4(7M(r>) = 4(7) E C ,, *^W^(»)lWx„W 

E ^ W '(7)1(^(»)1WX.» 



(13) 



The expansion coefficients j^ xn y nz > Sz 
of the basis vectors: 



(7) are calculated by using the one-dimensional rotation 

^HnM^nM = E BnlnK^Sx) Tpn y (v), (14) 



n y = n' x + n' y , are calculated 



and the coefficients Bn^nK'j), which are non-zero only for n 2 

numerically by using the Gauss-Hermite (G-H) quadratures. Rotations around the y— and z 
axes by the Euler angles j3 and a, respectively, are performed in exactly the same way. This 
method is exact, numerically very efficient, and inherent to the hfodd code. Indeed, since 
the rotated state is expanded in the original HO basis the structure of the code is preserved 
and allows, after only a minor generalization, for using the existing subroutines of the code to 
calculate the Hamiltonian Tl(Q) and norm jV(fi) kernels. 

Calculation of the kernel, (<&|0|$), for an arbitrary quantum operator O, proceeds along 
standard rules for calculating matrix elements between two non-orthogonal Slater determinants, 
see e.g. [16]. In particular, the norm kernel is given by: 



($|$) = Det O, 
with the overlap matrix elements defined as: 



Oi 



V 



dr 



,ifi(r,a)(p*(r,a). 



(15) 



(16) 



2 The present implementation assumes a spherical HO basis to ensure that the rotation does not induce 
components of single-particle wave functions that would go beyond the assumed basis. 
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The kernel of an arbitrary one-body operator F reads: 



= Y,(Vi\m)0^ = JJ dfdf'Y: (m|F|fV) p(r W,fa), (17) 

where OZ denotes the matrix element of the inverse of the overlap matrix O while p(f'a', fa) 
is the one-body transition density matrix defined as: 

p(f V, fa) = J2 ti(f, a) a') Oj. (18) 

ij 

At the same time, the most general GCM kernels can be calculated by using two different sets 
of wave functions for ipi(f,a) and (pj{f,a), which correspond to two arbitrary different Slater 
determinants. 

The kernel (TT7T) and density matrix ([18]) can be further simplified by introducing auxiliary 
ket-states defined as: 

^(rV)^^^,(fV)0^. (19) 
j 

Indeed, this allows for rewriting the transition density matrix to a "diagonal" form: 

p(f, a, f, a') = ]T <p*{f, a) a'), (20) 

i 

where the summation goes over a single index in full analogy to the diagonal HF density matrix. 
It means that the transition density matrix and one-body kernels can be calculated by using 
the standard subroutines of the code. All what needs to be done is a substitution of the HF 
ket-state by the auxiliary ket-state (TT9T) . 

Similar property also holds for the two-body-interaction kernel, which preserves the func- 
tional form of an energy density functional (EDF) derived for this interaction by averaging it 
over the Slater determinant. Again, all what needs to be done is a replacement of the density 
matrix by the transition density matrix. In our particular case, the Skyrme EDF can be ex- 
pressed by using six local isoscalar (t = 0) and six local isovector (t = 1) densities, including 

the particle p t , kinetic r t , spin s t , spin-kinetic T t , current j t , and spin-current J t densities and 
their derivatives, see Eqs. (1-5)— (1-7) . The Skyrme- interaction kernel preserves the functional 
form of the Skyrme EDF: 

(*\Vsk\$) _ H Sk (tt) = 




+ Cl[p Q ]f t + ■ M t + cp t ■ f t + Cfj t + c7% ■ (v x I) . 

(21) 

— * — * *~ * 

It depends on six local transition densities p t , ft, s t , Tt,j t , Jt, which are the counterparts of the 
six diagonal local densities mentioned above. 
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However, in the case of density- dependent interactions, like the Skyrme force, the kernel is 
not fully defined. The procedure must be augmented by a prescription concerning the treatment 
of the density-dependent term. In the hfodd code we implement the standard prescription, 
that is, in the primary coupling constants we replace the isoscalar density by its transition 
counterpart, i.e., Cf [po] — > Cf [po] and C t s [po] — > C* [po] (see discussion in Ref. [T7]). 

As already mentioned, implementation of the AMP requires a relatively minor recoding 
of the standard routines of the code. This concerns two generic matrices, D and L (1-39), 
which are used in the hfodd code to encode density matrices and their derivatives. In version 
(v2.38j), these matrices were generalized in the following natural way: 

D%, a = E (V A <Mr>)) (V^*(rV)), (22) 

i 

and 

Lm ' = ^ E (AHM(^) + A0,(faM(rV)), (23) 

where = (1, V) for ft = 0,1,2,3 and the indices q and q' denote the signs of spins a and 
a', respectively. In order to allow for non-hermitian, complex transition densities and their 
derivatives, the set of formulas expressing them in terms of the D and L matrices, see Eqs. (I- 
40)— (1-51), was generalized in the following way: 

• scalar densities 

P = D++ + D--, 

Ap = 2f + 2(L ++ + L— ), 
V-J = i(D+ + -D^--D+ + + D^-)+i(D+- + D^+-D+ 2 
-(Dfr-Dj+-D&- + D£-), 

• vector densities 



h = D+- + Dw + , (28a) 

h = i(D+--Dw + ), (29a) 

h = D++-Dm, (30a) 

f i = + ( 31a ) 

f 2 = ^K-V). ( 32a ) 

f 3 = -!>-)> ( 33a ) 



(24) 
(25) 

(26) 
(27) 
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Asi = 2Ti + 2 (L + - + , (34a 

As 2 = 2T 2 + 2i(L + ~ - L~ + ), (35a 

A§ 3 = 2T 3 + 2 (L ++ -L— ), (36a 

v m p = d++ + d;- + d++ + d o7 , (37 

^ = ^D; o ++D- --D++-D^-, (38 



(V = (d++ - D^T + D# " ^cf ) - <^oV + Dt - - Dj? - J D 3 o + ),(39a 

(V x s) 2 = + !>"+ + 2?+" + L>3 + ) - (l>++ - D Q T + L>+ + - ^-),(40a 

(V x s) 3 = i(D+- + D+- - D + - £>r+) " (^oV + + D 02 + + Ifc+^la 



(Vx4 = i{D^-D^ + D^--D^), (42a 

(Vxj) 2 = i (D++ - D+ + + D~ -L> lT ), (43a 

(Vxj) 3 = ^^V-^ + ^iT-^T), (44a 

• tensor density 

J* = ^{D^ + D^-D^-D-+), (45a 

^2 = l(Dfr-Drf-Dfr + D£), (46a 

^3 = ^{DU-D-.o-D^ + D,-). (47a 



The final step in the calculation of the Hamiltonian and norm kernels flU)-© is the numerical 
integration over the Euler angles. In our implementation, these three-dimensional integrals are 
calculated by using the Gauss quadratures [18]. To achieve a maximum accuracy, the Gauss- 
Tchebyschev (G-T) quadrature is used for the integration over the a and 7 Euler angles and the 
Gauss-Legendre (G-L) quadrature is used for the integration over the (3 angle. This combined 
technique ensures, in fact, exact integration, provided that the numbers of the G-T nodes 
n a = n 7 and the number of the G-L nodes rip are sufficiently large. 

The numbers of the Gauss nodes to be used depend on the value of the maximum-spin 
component I ma x m the Slater determinant expansion in good angular-momentum basis states 

1 max I 

\*)= E E \IKK). (48) 

I — 1-min K — / 
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Therefore, condition 2I max < rii — 1 for i = a,/?, 7 is common for all three Euler angles. In 
practice, twice as large numbers of nodes ensure good numerical stability for terms that are 
not polynomial, like for example the density-dependent terms. 

Note, however, that since I m ax is a priori unknown, the precision of integration can be 
verified only a posteriori. The user can verify this precision by inspecting the completeness 
relations for the overlap: 

1 = Y,(IKK\IKK) = J2 / dnD^imm, (49) 

IK IK ° n J 

and the Hamiltonian (or the Hartree-Fock energy Ehf)' 

E HF = ($|tf|$> = Y.( IKK \H\IKK) = E^irir / dMD% K ($\HR(n)\$). (50) 

IK IK ° n J 

The results concerning the completeness relations constitute part of the standard printout of 
the code. 



2.2 Calculation of electromagnetic transition probabilities 

Angular momentum projection opens up a possibility to compute fully quantum mechanically 
transition rates for electromagnetic radiation between final (f;IfMf\ and initial |i;/jMj) K- 
mixed states (j2J) of angular momenta If and I iy respectively. Reduced transition probabilities 
are defined by means of the reduced matrix elements, which are provided by the code, as : 

B(EX, U — > If) = -J—^lfllQ^l^, 

fl(MA,I,— ►!/) = ^-^I^J/IIMaIK;/,)! 2 , (51) 

for electric = r x Y\n(cp, 9) and magnetic = g s M s -\^ + 9iMi;\n transition operators, 
where g s and gi denote spin and orbital gyromagnetic factors, respectively. The spin and orbital 

parts of the magnetic transition operator equal M s; a m = (VQa^) ' & an d = ^-(VQa/J • L, 
respectively. 

The electric and magnetic transition operators transform under spatial rotations 
like components of spherical tensor of rank A. According to the Wigner-Eckart theorem, 
the matrix elements of these operators are, therefore, equal to: 



(/;/ / M / |T AM |z;/ i M,)- v -^ ^ hM ^ 

fZlf + l 



« (hlgmM, (52) 



where Cjj M . f x denotes the Clebsch-Gordan coefficient. 

The left-hand side of Eq. (152"]) can be calculated by using the projected states (pQ). Using 
the definition of the i^-mixed states (TSl), it can be written as: 



(f;I f Mf\f x ^;IM) = E 9 { I^hK^f M f K f\ f ^ M ^) 

KtKf 

E 9^ f 9^f\P^M f fx,Pk Ki m, (53) 



KiK 
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where \<$>f) and are, in principle, different Slater determinants. The matrix element (153!) 
can be further simplified after applying the following transformation rule, valid for arbitrary 
spherical tensor: 

p 7 / f, pU _ n I f M f \ - / 1 ^2n' r< I f K t T P 7 * (KA\ 
r K f M f 1 ^ r M i K i — ° hMi \n V L ) ^hM Xjj,' 1 MK,- \°^) 

Indeed, by using this formula one obtains: 

(ifMfKfif^iMKt) = egg, EC- 1 )^'^^^/^^!**). ( 55 ) 

what brings the calculation of the transition rates down to a standard task of the AMP, namely 
to the calculation of the matrix elements of the one-body operator {^f\T\fjiP^ K . which is 
described in detail in Section [27T1 By comparing Eq. (1521 for / = fT/ and % = Ki with Eq. fl55|) . 
we finally have: 

(/^HfkHJ^) = (•f/lf^P^J^). (56) 

These quantities enter directly the equation for the reduced matrix elements between the AMP 
X-mixed states, which reads: 

(/jJ/HTkH*;/*) = J2 iiX&ilfKfllTMKi). (57) 

K t K, 

For the sake of completeness, it should be recalled that in the code hfodd the multipole 
operators are defined as Qx^ = r x Y Xfl , see Eq. (IV-2). Hence, the matrix elements of multi- 
pole electric and magnetic operators are calculated in the code for dual tensors. In order to 
be consistent with the matrix elements of multipole operators printed by the code one needs, 
therefore, to take T\ M = Q* x for electric transitions and use VQ* Xfl in the definition of mag- 
netic transitions. Of course, these phase conventions do not influence the reduced transition 
probabilities (EH). 

2.3 The Yukawa interaction 

The code hfodd can now evaluate the expectation value of a two-body Yukawa interaction 

g~ ar ij 

V = 22 [a + hoi ■ <jj + CTi ■ tj + doi ■ a k n ■ r k ] , 

i<j r ij 

where = |r$ — r ^ | , and a, a,b,c, and d are arbitrary constants. The code makes use of an 
approximate expansion of the Yukawa function in terms of Gaussians: 

— « 6.79 e" 34x + 2.41 e ~^ x \ 0.786 e~ L44x + 0.241 e -°- 38:c - 0.062 e" - 15 * + 0.078 e' - 13 ^ , (58) 
x 

Of course, the Yukawa function is singular at the origin and Gaussians are not, but the volume 
element contains a factor of r 2 that makes the approximation quite accurate near the origin (it 
is a bit less so at larger distances). 
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As described in Ref. [19], hfodd can also use a series of Gaussians to approximate the 
expectation value of the CP-violating potential produced by pion exchange with a CP-odd 
pion-nucleon vertex: 



g rn^h c 



Vpr(ri - r 2 ) = — — - — { (tr 1 - er 2 ) ■ (r x - r 2 ) g Q n ■ f 2 - — (t 1z + t 2z ) 



9i 



2 

+92{3t 1z t 2z - fi ■ f 2 ) - y (o-i + cr 2 ) • (ri - r 2 ) {r u - r 2z 



exp(-m n \r 1 - r 2 \) 
x — : — — 



m^ n - r 2 
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1 + 



m 7T \r 1 - r 2 \ 



(59) 



where the g's label the isoscalar, isovector, and isotensor pion-nucleon coupling strengths and 
m-x and tjin stand, respectively, for the pion and nucleon mass in units of fin -1 . To simulate the 
effects of short-range correlations, absent from HF wave functions, we include a phenomeno- 
logical correlation function [20] 

/(r) = 1 - e~ Llr2 (1 - 0.68 r 2 ) , (60) 

and make the approximation 
g (r) = f( r f —— (! + -!_)« 1.75e" Llr2 +0.53 e - - 68r + 0.11 e~°- 21r \ 0.004 e" 06 '' 2 , (61) 



r 2 \ m^r 

where m n = 0.7045 fm _1 and the numbers in the fit all have units of fm -2 . The extra factor of 
r not included in Eq. (IBTj) (i.e. the factor r\ — r 2 in Eq. (159^1 ) is treated separately. 

A Gaussian interaction between, e.g., particles 1 and 2, leaving out the spin dependence, 
factors into functions that each depend only on the x, y, or z component of the vector difference 
ri — r 2 . To evaluate the expectation value of the first Gaussian, which depends on x\ — x 2 , it is 
easiest to write the product of two oscillator wave functions (of x% and x 2 ) in terms of relative 
and center-of-mass coordinates. To do so, we use 

f _ Q l + °4 f _ Q l ~ a 2 (ao \ 

a CM — ~y| ) a rel — ? l DZ J 

to obtain the relatively compact one-dimensional "Moshinksky bracket" 

^n\N\ ni \n 2 \ (_l)n-n 1+ i 

(ni, n 2 \n, N) = b nl+n2>n+N —— — — — - , (63) 

2 2 i i>-{N — l )-\ n i — iy\n — ni + 1)\ 

where n\, n 2 refer to one- dimensional oscillator wave functions associated with particles 1 and 
2, and n, N label oscillator wave functions in the corresponding relative and center-of-mass 
coordinates. The expectation value of the interaction then involves the integrals of Gaussians 
of arbitrary range against two Hermite polynomials, with additional factors of x, y, or z if the 
T-violating potential in Eq. (1591) is being used. The resulting expressions are straightforward. 

Once the integrals have been evaluated for each Cartesian coordinate, they must be folded 
with the two density matrices to give the final expectation value. To do this, the code must 
sum over a large number of individual oscillator quantum numbers. Explicitly (and in general), 
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the expectation value of a term in the Gaussian approximation to a finite-range potential like 
the ones we have been considering reads 
rpr'k'rk \ " \ " V r r<x rry nz n k! (m' n' \rr k (m n \ 

-^dir l—i r 'nJ 's',m 'W ' h ns ,mu K -' m' x m x ,n' x n x m y m y ,n y n y ^ m' z m z ,n' z n z u u' s' \' ' b y' b yJ u us\' ' b y' ''V J • 

m'n'mn s'u'su 

(64) 

Here p T ns mu is the density matrix for particles of type r=p,n in the basis specified by the 
oscillator quantum numbers n=(n x n y n z ) or m=(m x m y m z ) in the three Cartesian directions 
and the simplexes s or u, G" 1 , , are one-dimensional integrals representing matrix elements 

of a Gaussian (or a Gaussian times a coordinate) in the i=x, y, or z directions, and &ts( m y n y) 
for k=0, x, y, or z are the matrix elements of the Pauli matrices (with a = 1) in the simplex 
basis (1-78) that depend on the numbers of quanta in the y direction (m y n y ). 

Altogether, sums in Eq. (1641) require 12 independent sums over the oscillator states, i.e., 
about A^ 2 operations for a basis cut at A" oscillator shells, and are intractable unless performed 
efficiently To reduce the number of operations, hfodd sums over the simplex quantum numbers 
first, resulting in the intermediate quantities 

-^raW = X! Pn>s',m'u>< J t>s>( m y n y)> (65) 
s'u' 

D nm = P™,mu°Zs( m vnv) > ( 66 ) 

su 

the computation of which requires A"q operations. After these are stored, the sums in each 
Cartesian direction are done one-by-one. hfodd first computes auxiliary ^-direction matrices 

yrk _ \ " rjrfe QV (F>7) 

n x m' y n z ,m x n' y m z / , n x n y n z ,m x m y nn z m' y m y ,n' y n y > V u / 

a task that requires A"q operations. Next it performs the z sums, yielding 

n x m y m' z ,m x n' y n' z / j 1 n x m l y n z ,m, x n' y m z r m' z m z ,n' z n z 1 V uo / 
n z m z 

the computation of which again requires A"q operations. The sum in the x direction is postponed 
to save storage (the auxiliary matrices Yn x m' n Z) m x n' m z an d ^ntm' m' ,m x ri ri are eacn f° ur ti mes 
larger than the original density matrix). Instead, for fixed values of the x quantum numbers n x , 
n' x , m x , and m' x , the code first sums over the "primed" y and z labels, computing and storing 

v-r'k',rk _ \p n T ' fc ' 7 rk 

n' x n x ,m' x m x / t n^n',ni,m^mi,mi n x mi,rnL,m x ni,Ti', ' v / 



-y-z 



a procedure that once more requires A"q operations. The final result is obtained by summing 
over the Nq x-direction labels, giving 



pr'k'rk _ y-7 -'A ,rk (~<x (7C\\ 
dir / j n' x n x ,m' x m, x (j m l x rn x ,n' x n x " 

n x n x ,m' x m x 

The separation of the Gaussian interaction into its Cartesian components is what allows the 
original A^ 2 terms to be reduced to Nq terms. For comparison, the mean-field matrix elements 
calculated for a zero-range interaction require a sum over about Nq terms, see I. 

Figures [T] and [2] show the actual CPU times required in calculations using the Skyrme 
and Yukawa interactions, respectively. At Nq = 10, the latter take more than two orders of 
magnitude longer and, moreover, their CPU times scale like Nq rather than A"q . A local EDF 
clearly makes for easier computing. 
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Figure 1: The hfodd CPU times required for calculations that use the standard Skyrme EDF, 
shown as a function of the number of HO shells N . The doubly logarithmic scale in the Figure, 
shows that these times scale as N%. 




N 

Figure 2: Same as in Figured] but for calculations that use the Yukawa interaction. The CPU 
times scale as NL 
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2.4 The BCS solutions for state-dependent pairing gaps 

In the previous versions of the code hfodd, only the seniority pairing force has been imple- 
mented in the particle-particle channel. In version (v2.38j), we have introduced the subroutine 
delpai which solves the BCS equations for the seniority pairing force, the seniority pairing 
force with fixed pairing gap, as well as for the zero-range density-dependent pairing force. The 
latter option amounts to using the state-dependent pairing gaps. 

For the seniority pairing force, the neutron and proton pairing strength constants are defined 
as in Ref . [21] . These values can be scaled by the multiplicative factors FACTGN and FACTGP for 
neutrons and protons, respectively (see Section II-3.2). 

The zero-range density- dependent pairing force, Eq. (IV-14), is defined by the sets of param- 
eters {Vq, V\, a} or equivalently by {Vo, po, a}, see Section IV-3.1. The pairing matrix elements 
between pairs of time-reversed states (the array GMATRI) are calculated in subroutine ginter. 
After solving the BCS equations for state-dependent pairing gaps A k , the code calculates the 
average neutron and proton (spectral) gaps defined as: 

were v k and Uk are the BCS occupation amplitudes. 

2.5 Calculation of Bohr deformation parameters 

The shape of the nuclear surface is usually described in terms of the multipole expansion, 

oo A 



R(6, <f>) = R 1 + £ £ axM9, 0) . (72) 

y a=ia*=-a / 

In self-consistent methods, the radius R and the deformation coefficients a\^, A > 1, have to 
be determined on the basis of the calculated density distribution according to some convention. 
The general concept is to consider a sharp-edge body of constant density po and shape defined 
by Eq. (jT2"l) . and choose po, Ro, and so that the monopole surface moment, Qq (IV-3), and 
the electric multipole moments, (IV-2), A > 0, of the body equal those of the mean-field 
solution. 

The moments Qq and of any density distribution, p(r), are defined as average values 
of the functions given by Eqs. (IV-2) and (IV-3), i.e., 

Q s 00 = J d 3 fpr 2 , Q x , = a^J d 3 f pr x . (73) 



The conventional factors are listed in Table IV-5; in particular, aoo = v47r. The monopole 
components, Qq and Qoo, are related to the root-mean-square radius, R rm s, Qoo = Qoo-R? ms , 
and particle number N, Qoo = N. Below, the quantities pertaining to the mean-field solution 
are marked with bars, R Tm s, Qoo an d Qx^- For the reference body, Eqs. (1731 take the form 

Q s 00 (po,R ,a) = ^ I** &<P (\medeR\e )( p) , (74) 
5 Jo Jo 

Qx^Po,Ro,a) = a AM ^| o 27r d0^sin0d^ A+3 (^,0)V A ;(^,0) , (75) 
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where a denotes the ensemble of the parameters «a m - Thus, the equations of the problem read 

Qo (po, Rq, a) = , Qa m (po, Ro, a) = Qa m • (76) 

They cannot be solved analytically. 

For small deformations, however, one can determine po and Rq from Qq and Qoo by setting 
c^X/j, = in Eqs. (1731) and (175|) . Then, one can find by expanding the expression in Eq. (1751) 
up to first order in about ct\^ = 0. This leads to the linear approximation, 

3 Qoo [5 - 4vr Qa m fvv , 

P °-^lf' R ° = i3 R ™> ^-S^Qo^' (7?) 

which was used in previous versions of the code hfodd, see Section III-2.9. 

In version (v2.38j), exact solutions to Eqs. ( 1761) are sought numerically, for arbitrary defor- 
mations. First, the integrals in Eqs. (174|) and (1751) are evaluated by using the G-L quadratures. 
Then, the problem is formulated in terms of fixed-point equations, which are solved by using 
the standard iterative method. Iterations stop when Eqs. (176|) are satisfied up to an accuracy 
which guarantees that all the digits printed on the output are correct. It may happen for 
very large deformations that the procedure diverges and the exact values of the deformation 
parameters are not found. 

Unlike in the linear approximation, the values referred to as exact depend on the value of 
the maximum multipolarity A max considered in Eqs. ( 1761) . This is not a very dramatic effect, 
and in order to estimate the exact values for a multipolarity A, it is recommended to perform 
calculations for multipolarities up to the next one of the same parity, i.e., for X max = A + 2. 
In the code, X max is set independently of the number of multipole moments calculated for the 
mean-field solution, although it must not exceed the latter one, of course. 

According to the adopted convention, see Section III-2.9, the code prints the Bohr defor- 
mation parameters, /?a m , which contain an additional factor of \pl introduced for p ^ 0, 



^ = p - <V) • (78) 

2.6 Constraints on the Schiff moments and scalar multipole moments 

Apart from calculating and constraining the multipole moments Qa m (IV-2) and surface mo- 
ments Qf M (IV-3), the code hfodd version (v2.38j) also calculates the average values of the 
so-called Schiff multipole moments, 



Ql(r)-Ur 2 )Qx»(r), (79) 



Ql(r) = %r-f(rV 

where for the neutron, proton, or total Schiff moment, (r 2 ) denotes the neutron, proton, or 
total mean-square radius, respectively. Note that the Schiff moments depend through (r 2 ) on 
the self-consistent solution; therefore, in the given iteration of the code, values of (r 2 ) are taken 
from the previous iteration of the code. In this way, a proper definition of the Schiff moment is 
ensured only after the convergence is reached. Note also that the factor of ^, which is included 
in the standard definition of the Schiff moment [221 [19] is not included in the definition of 
Eq. Q7J). 
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In version (v2.38j), the code hfodd can calculate and constrain either the surface moments 
or the Schiff moments, according to the switch ISCHIF, see Section 13.51 below, but not both 
simultaneously. 

Since the physical constraints on the Schiff moments have to be set for the proton moments 
only, the meaning of parameters IFLAGS (see keyword SURFCONSTR in Section IV-3.5) has been 
generalized in the following way: As previously, for IFLAGS=1, the constraints on the surface 
or Schiff moments pertain to the total moments, but for IFLAGS=2 or 3 they now may pertain 
to the neutron or proton moments, respectively. 

In version (v2.38j), the code hfodd can also set constraints on scalars Q\ built from the 
multipole moments Q\^. 



Qx — oao. 



E (80) 




For A = 2, such a constraint is useful when the minimum of energy is to be found as function 
of the deformation 7 for a fixed value of the deformation (3. Constraints on scalars are also 
insensitive to the overall orientation of nucleus in space. 

2.7 Quasiparticle blocking for the HFB solutions in odd and odd-odd nuclei 

In order to include pairing correlations, hfodd solves the HFB equations: 

w ( x <p) = (x v) 

where 

"-(-* -A a) m 

is the HFB matrix. Here, h' is the s.p. Routhian operator, A is the Fermi energy, A is the 
antisymmetric pairing potential, and E is a diagonal matrix of quasiparticle energies. In the 
case of even-even nuclei, \ (resp. f) are the wavefunctions of quasiparticles with positive (resp. 
negative) energies. Note that these are 2M x M matrices, where M is the dimension of the s.p. 
basis. One may write them down in terms of upper and lower components using the M x M 
matrices A and B: 



*={b)' *={a>)- < 83) 

The matrices A and B form the matrix of the Bogolyubov transformation: 



In the case of the so-called proper Bogolyubov transformation, for which det^4 = 1, the HFB 
state is a mixture of states with even number of particles only. The blocking of quasiparticle 
k is done by replacing one column in the tp matrix by the column which corresponds to the 
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quasiparticle state with the opposite energy from the \ matrix. This leads to exchange of the 
two columns of the Bogolyubov transformation matrix A and det.4. = — 1. Now A corresponds 
to the improper Bogulyubov transformation and the resulting HFB wave function is a super- 
position of states with odd number of particles only [H]. Let Ai denote the i-th column of the 
matrix A. Then the corresponding HFB wave functions of blocked states have the form: 
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(85) 
(86) 



To decide which quasiparticle should be blocked the code calculates the overlap of s.p. and 
quasiparticle states in the following way: overlaps of the s.p. state with the upper component 
and with the time-reversed (complex-conjugate) lower component of the quasiparticle state. 
The greater of these two is chosen. The quasiparticle state with the largest overlap with the 
selected s.p. state is chosen to be blocked during each iteration. 

In the case of conserved simplex symmetry, the Routhian and pairing potential have the 
form: 



ti 




A 




17) 



Then, the general HFB equations decouple into two different set of independent equations. The 
code hfodd solves the equation: 



ti + — A A+ 
-A* til + A 







-E 4 



from the solution of which, one may reconstruct the complete solution of Eq. (18T|) : 
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(89) 



Blocking of the k-th quasiparticle (for k < M/2) leads to the following expressions for the 
wave HFB functions (|85|) and (|86l): 
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where horizontal and vertical lines are introduced to distinguish M x M/2 blocks. One can 
easily derive analogous expressions in the case of M/2 < k < M. Instead of using the blocked 
matrix ip in the form of (1911) . the code uses instead a 2M x (M + 1) matrix ip': 
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(92) 



It can be easily seen that this matrix leads to the same generalized density matrix 1Z = <fi(p'. In 
this representation, the wave function of the blocked quasiparticle is just zeroed and the wave 
function of the quasiparticle with opposite energy is added as an additional column. 

2.8 The Broyden method 

The HFB equations are a system of non-linear integro-differential equations, which require to 
be solved self-consistently: an initial guess for the density matrix and pairing density is used 
to generate the HF potential and pairing field and define the HFB matrix; Solving the HFB 
equations yield a new set of eigen-functions that are used to calculate the densities at the next 
step, and this loop is executed until convergence is achieved. Different criteria can be used for 
convergence. In hfodd, iterations are stopped when the difference between the HFB energy 
and the sum of s.p. energies is less than EPSITE, see I and II. 

Mathematically, the self-consistent HFB equations can be viewed as a particular example 
of a fixed-point problem. Formally, they read: 

Vffl = I(vL m) ), (93) 

where V;™' 1 is a vector of size N containing certain initial conditions (i.e., the set of wave- 
functions, HF potential, matrix elements of the HFB matrix, etc.) at iteration number m. The 
solution V to the HFB equations satisfies: V = I(V), or: F(V) = V - I(V) = 0. By default, 
most of nuclear structure codes iterate the vector V by taking as input to iteration m + 1 a 
linear combination: 

vL w+1) =aVS + (l-a)vH. (94) 

This so-called linear mixing scheme was implemented in previous versions of hfodd. 

If the function F(V) is differentiable, the roots of the equation F(V) = can be found by 
a generalized Newton-Raphson method: this is the basis for the Broyden method to accelerate 
convergence. At each step m of the self-consistent process, the input to the next step is 
computed from: 

V £H- 1) =V M_ B (m) F (m) j (95) 

where B^ m ^ is an approximation to the inverse of the Jacobian of the vector function F(V) at 
iteration m. The exact expression of Eq. (1951) was first given in [23]. The modified Broyden 
method is a slight modification of the original Broyden method designed to avoid storing (po- 
tentially large) N x N matrices. It was presented in [21] and applied in a variety of nuclear 
structure problems in [25]. In practice, it is a correction term added to Eq. (1941) . 
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The modified Broyden method was implemented in hfodd according to the formulation of 
[25J. In hfodd the quantities that are iterated are the components of the HF fields on the G-H 
integration mesh, see I. The Broyden vector therefore contains the value of the components of 
the self-consistent HF fields at the G-H node (xi, yj, 2%). The size N of this vector is proportional 
to the total number of field components and the numbers of G-H nodes along each direction 
NXHERM, NYHERM , NZHERM. It is given by: N =44x(NXHERM x NYHERM x NZHERM). For the 
Broyden method applied within the Lipkin-Nogami (LN) calculations, the matrix elements of 
the density matrices have to be stored in addition, which adds 2 x M 2 elements to the Broyden 
vector (M being the size of the s.p. basis). 



2.9 The Lipkin-Nogami method 

Methods of restoring the particle-number symmetry must be implemented in studies of pairing 
correlations because some observables like even-odd mass staggering or pair-transfer amplitudes 
are influenced significantly. In addition, the quantitative impact of the particle-number restora- 
tion depends on whether the pairing correlations are strong (open-shell systems) or weak (near 
closed shells). 

The LN regime of hfodd gives an efficient way of performing approximate particle-number 
projection (PNP) calculations. It is based on the LN method [261 considered as a variant of 
the second-order Kamlah expansion [2"5j |2"9"| I3"0"| I3"T] , in which the PNP energy is approximated 
by a simple expression, 

£tOT = &HFB + ^LN, (96) 

where Shfb is the HFB energy, and £ln is the LN correction, 

£ln = -A 2 ((iV 2 ) - N 2 ) = -2A 2 Trp(l - p), (97) 

with A 2 depending on the HFB state and representing the curvature of the energy with respect 
to the particle number. In the Kamlah method A2 is varied along with variations of the HFB 
state while in the LN method this variation is neglected and A2 is simply evaluated after each 
iteration in order to find the best estimate of the energy curvature. 

When the HFB method is applied to a given Hamiltonian, values of A2 can be estimated 
by calculating modified mean-field potentials that are analogous to the standard HFB mean 
fields. However, in the spirit of the EDF approach, in hfodd is adopted an efficient algorithm 
[32] estimating the curvature A2 from the seniority-pairing expression, 

_ G cS Tr'(l - p) K Tr> - 2 Tr(l - pfp 2 
2 4 [Trp(l - p)f - 2 Trp 2 (l - p) 2 ' 1 ' 

where Tr„4 = J2 n A nn stands for the trace of a matrix A, while Tr'^4 = J2n-^nn- The effective 
pairing strength, 

A 2 

G cff =--^, (99) 
is determined from the HFB pairing energy, 

£ P air = -^TrA K *, (100) 
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and the average pairing gap, 

* Tr'Ap , , 

A = ^r- < 101 > 

All quantities defining energy £ LN , Eq. (157)1 . and A 2 , Eqs. (1951)- (j 101ft . depend on the self- 
consistent solution, and the microscopic interaction. In hfodd, they are calculated in the 
subroutine LIPCOR using the canonical representation for the density matrices in terms of the 
canonical occupation amplitudes L4 and 14- 

As a result of the LN calculations, pairing correlations never collapse, which is also the case 
of the exact PNP before variation. A comparison of different HFB particle-number projection 
results is further discussed in Ref. 132 1. 



2.10 Exact Coulomb exchange 

In evaluating the exact Coulomb exchange mean fields and energies we use the methods de- 
veloped for the Gogny interaction, as described in Appendix A. 5 of Ref. [33]. They are based 
on expanding the Coulomb interaction into a sum of Gaussians, which allows using the same 
infrastructure as developed for the Yukawa interaction discussed in Section 12.31 In particular, 
we use the identities: 

where the second integral has been reduced to a finite domain by the substitution a = 
£ (1 — £ 2 ) and L stands for the largest of the three HO lengths L M = ^Jh/mu^, \i = x,y, z. 
By using the G-L quadratures we can now present the Coulomb interaction as a finite sum of 
Nc Gaussians: 

- = Y^A i exp(-a i r 2 ), (103) 
r i=i 

where constants A{ and a, depend on the G-L weights Wi and nodes as 

A ^B^-^ 3 ' 2 • a ^TH^w> (104) 

It turns out that this method provides for very precise and rapidly converging results for 
the exact Coulomb energies. This is illustrated in Fig. [31 where the error in the exact Coulomb 
exchange energy, AE exc , is plotted as function of the number of Gaussians Nc- A quite precise 
estimate is obtained for Nc = 7 Gaussians and the machine precision is obtained by doubling 
this number (along with doubling the CPU time). For Nc = 7, these CPU times are shown in 
Fig. HI which illustrates the fact that these times scale with the number of HO shells as iVj. 

2.11 Utility options 

The following utility options have been implemented: 

1. In all the keywords, the minus character "-" can always be used in place of the underscore 
character "_" and vice versa. 
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Figure 3: Errors in the exact Coulomb exchange energy, AE CXC , plotted in function of the 
number of Gaussians N c used in Eq. (I103I) . 
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Figure 4: Same as in Figure [T] but for the calculations that determine the exact Coulomb 
exchange energies. The CPU times scale as Nq. 
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2. An external HO potential can be added to the self-consistent mean field to calculate 
properties of an atomic gas in a trap, see the keyword INSERT_H0 below. 

3. Arbitrary values can be added to the coupling constants, see the keywords EVE_ADD_TS, 
0DD_ADD_TS, EVE_ADD_PM, and 0DD_ADD_PM below. 

4. Calculations with fixed values of the Fermi energies can be performed, see the keywords 
FIXFERMI_N and FIXFERMI_P below. 

5. For the HFB calculations, eigenvalues of the HFB mean-field s.p. Hamiltonian or Routhian 
can be printed, see the keyword HFBMEANFLD below. 

6. Within the HF calculations, the filling approximation can be used, see the keywords 
FILSIGJJEU and FILSIG_PRO below. 

7. Constraints on the intrinsic spin only can be used, see the keyword NORBCONSTR below. 

8. Various one-line data can be printed during the iteration, see the keyword 0NE_LINE 
below. 

9. Integrals of several symmetry-violating terms can be printed, see the keyword PRINT_VIOL 
below. 

10. The Nilsson labels with respect to the x, y, or z axes can be printed, see the keyword 
NILSSONLAB below. 

11. Matrix elements of the mean-field Hamiltonian can be saved on the disk, see the keywords 
FIELD_SAVE, FIELD_0LD, REC_FIELDS, and CONTFIELDS below. 

2.12 Corrected errors 



2.12.1 Euler angles. The Euler angles have been calculated for the quadrupole moment, i.e., 
for the complex conjugate spherical harmonics, instead of the complex conjugate quadrupole 
moment, i.e., for the spherical harmonics themselves. 

2.12.2 Predefinition of STIFFA. In subroutine PREDEF, the value of variable STIFFA has been 
incorrectly predefined to 0.01 instead of 0.00. 

3 Input data file 

The structure of the input data file has been described in Section II-3; in the present version 
(v2.38j) of the code hfodd this structure is exactly the same. All previous items of the input 
data file remain valid, and several new items are added, as described in Sees. I3.1H3.8I 

Together with the FORTRAN source code in the file hf odd.f , several examples of the input 
data files are provided. File dyl52-f . dat contains all the valid input items, and the input data 
are identical to the default values. Therefore, the results of running the code with the input 
data file dyl52-f.dat are identical to those obtained for the input data file containing only 
one line with the keyword EXECUTE. 
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File ge064-a.dat contains the input data that allow for determining the rotating triaxial 
state in the nucleus 64 Ge and then performing a 3D AMP of this state. This data file comprises 
four runs of the code: (i) starting from the Nilsson initial potential, the code performs 20 
iterations with constrained values of quadrupole moments Q20 and Q22, such that the initial 
deformation and orientation of the nucleus is determined, (ii) after releasing the constraints, the 
triaxial minimum is converged, (iii) the rotating solution for / = 6 is converged by setting the 
angular frequency of tkv = 0.575 MeV, and (iv) the AMP of the rotating solution is performed. 
Small numbers of G-T (10) and G-L nodes (10), see Section |2TT| which are used in this run, do 
not allow for a precise AMP resolution of high angular momenta. This data file is only meant 
to provide an example of rapid calculation. 

File ge064-b . dat contains the input data that allow for a correct AMP, with higher numbers 
of G-T (40) and G-L nodes (40), but its execution requires a CPU time which is 4 3 =64 times 
longer (about 200h). Alternatively, one can run in parallel 40 jobs by executing the input data 
files ge064-c.dat with characters NUASTA replaced by integers from 1 to 40. 

File snl20-b.dat contains the input data that allow for determining the ground state in 
the nucleus 120 Sn with the Lipkin-Nogami corrections taken into account. 

File ge064-a.dat is reproduced in section TEST RUN INPUT below. Files ge064-a.out 
and snl20-b.out contain results of executing code hfodd version (v2.38j) for the two corre- 
sponding input files. Selected lines from the output file ge064-a . out are reproduced in section 
TEST RUN OUTPUT below. 

3.1 Interaction 
Keyword: C0UL0MBPAR 

7, 1, 1 = IC0TYP, IC0UDI, IC0UEX 
For IC0UDI=0, 1, or 2, the Coulomb direct energy and Coulomb mean field are neglected, 
calculated by using the Green- function method, see Section 1-5, or calculated by using the 
Gaussian-expansion method, see Section I2.10[ respectively. Similarly, for IC0UEX=0, 1, or 2, 
the Coulomb exchange energy and Coulomb mean field are neglected, calculated by using the 
Slater approximation (1-19) or calculated exactly by using the Gaussian-expansion method, 
see Section 12.101 respectively. For the Gaussian-expansion method, that is, for IC0UDI=2 or 
IC0UEX=2, positive values of IC0TYP denote the numbers of G-L nodes used in the integral 
of Eq. fU03p . For IC0UDI=2 or IC0UEX=2, the iteration can later be smoothly continued 
(IFC0NT=1, see the keyword C0NTFIELDS) only by saving the matrix elements of the mean 
field, that is, by requesting IWRIFI=1, see the keyword FIELD_SAVE. Therefore, IC0UDI=2 or 
IC0UEX=2 and IC0NTI=1 requires IFC0NT=1. 

Keyword: PAIR_MATRI 

1, 0, 0, = IDESTA, IDEMID, IDEST0, IDEDIS 
For IDESTA=1, the pairing matrix elements, required for the BCS pairing calculations with 
state-dependent pairing gaps (IPABCS=3), are calculated in the first iteration. At present, only 
the value of IDESTA=1 is allowed, because the option of storing the pairing matrix elements is 
not yet implemented. For IDEMID=1, the pairing matrix elements are calculated in the middle 
iterations, for IDEST0=1 in the last iteration, and/or for IDEDIS=n, in every n-th iteration. 

Keyword: INSERT_H0 
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= IPOTHO 

For IP0TH0=1, an external HO potential is added to the self-consistent mean field. Parameters 
of the potential are identical to those defining the HO basis. 

Keyword: EVE_ADD_TS 

0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. 

ARHCLT, ARH0_S,ARH0DT, ARH0DS, ALPR_T, ALPR_S, 

ATAU_T, ATAU_S, 
ASCU_T, ASCU_S, 
ADIV_T, ADIV_S 

By using this item, the coupling constants corresponding to a given Skyrme parameter set can 
be shifted by arbitrary values. The time-even coupling constants in the total-sum representation 
(1-14) are modified by adding the 12 numbers from ARHCLT to ADIV_S. By setting the scaling 
factors SRHCLT to SDIV_S equal to zero, see the keyword EVE_SCA_TS, one can input here a new 
set of the coupling constants. The name convention of variables is here the same as for many 
other variables in the code hfodd, see the keyword EVE_SCA_TS. 

Keyword: 0DD_ADD_TS 

0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. 

ASPI_T, ASPI_S,ASPIDT, ASPIDS, ALPS_T, ALPS_S, 

ACUR_T, ACUR_S, 
AKIS_T, AKIS_S, 
AR0T_T, AR0T_S 

Same as above but for the time-odd coupling constants. 
Keyword: EVE_ADD_PM 

0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. 

ARHCLP, ARH0_M,ARH0DP, ARH0DM, ALPR_P, ALPR_M, 

ATAU_P, ATAU_M, 
ASCU_P, ASCU_M, 
ADIV_P, ADIV_M 

Same as above but for the time-even coupling constants in the isoscalar-isovector representation 
(1-15). The total-sum additive factors are used first, and the isoscalar-isovector additive factors 
are used afterwards. 

Keyword: 0DD_ADD_PM 

0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. 

ASPI_P, ASPI_M,ASPIDP, ASPIDM, ALPS_P, ALPS_M, 

ACUR_P, ACUR_M, 
AKIS_P, AKIS_M, 
AR0T_P, AR0T_M 

Same as above but for the time-odd coupling constants in the isoscalar-isovector representation. 
Keyword: YUKAWATERM 

0.7045, 4.7565, 1.0, 0.0, 0.0, 0.0, 1, 

PIMASS, PNMASS, YUKAGT, YUKAGO, YUKAG1, YUKAG2, IYUTYP, I_YUKA 
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For I_YUKA>0, the code calculates the average values of the time-reversal- and parity-violating 
Yukawa interaction (|59|) . with the pion mass (m n ) of PIMASS and the nucleon mass (rajv) of 
PNMASS. If values of zero are read, variables PIMASS and PNMASS remain unchanged. Variables 
YUKAGT, YUKAGO, YUKAG1, and YUKAG2 correspond, respectively, to the coupling constants g, g , 
g~i, and <j2- For I_YUKA=2 or 3, the direct matrix elements of the Yukawa interaction (|59|) are, in 
addition, added to the self-consistent mean field. For I_YUKA=2 or 4, the exchange matrix are 
added. For IYUTYP=1, expression (16"TT) is used, while for IYUTYP=2, an analogous six-Gaussian 
expression is used without the short-range correction (I6UI) . that is, for f(r) = 1. For I_YUKA=0, 
all these input data are ignored and the Yukawa interaction is not taken into account. 

3.2 Symmetries 

Keyword: FILSIG_NEU 

2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, = 

KPFILG(0,0,0), KPFILG (0,1,0), KPFILG(1 , , 0) , KPFILG(1 , 1 , 0) , 
KHFILG (0,0,0), KHFILG(0,1,0), KHFILGQ ,0,0), KHFILG(1 , 1,0), 
K0FILG(0,0,0), K0FILG(0,1,0), K0FILG(1 ,0,0), K0FILG(1 , 1 ,0) 

These parameters govern the filling approximation for the neutron s.p. parity-signature con- 
figurations. Matrices KPFILG contain the indices of particle states in the four parity-signature 
blocks denoted by (+,+), (+,—), (— ,+), and (—,—), of given (paritysignature) combinations, 
i.e., (vr,r) = (+l,+z), (+1,— i), (— 1,+z), and (— 1,— i), respectively. Matrices KHFILG con- 
tain analogous indices of hole states, and matrices K0FILG contain numbers of particles put 
into the states between KHFILG and KPFILG, by using for them partial occupation factors of 
K0FILG/(KPFILG-KHFILG+l)/2. For K0FILG=0, in the given parity-signature block the filling 
approximation is inactive. The filling approximation is incompatible with pairing correlations, 
IPAIRI=1. 

Keyword: FILSIG_PR0 

2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, = 

KPFILG(0,0,1), KPFILG(0,1,1), KPFILG(1 , , 1) , KPFILG(1 , 1 , 1) , 
KHFILG(0,0,1), KHFILG(0,1,1), KHFILG(1 ,0, 1) , KHFILGQ , 1,1), 
K0FILG(0,0,1), K0FILG(0,1,1), K0FILG(1 , , 1) , K0FILG(1 , 1 , 1) 

Same as above but for the proton s.p. parity-signature configurations. 
Keyword: BCS 

-1 = IPABCS 

Parameter IPABCS defines the type of BCS paring calculations. For IPABCS=0, 1, 2, or 3, the 
BCS pairing calculations are, respectively, not performed, performed with the seniority pairing 
force, performed with fixed pairing gaps, or performed with the state-dependent pairing gaps. 
Value of IPABCS=— 1 is allowed for the sake of compatibility with earlier versions of the code, 
before (v2.13f). Then, IPABCS=0 is set for IPAHFB=1 and IPABCS=IPAIRI is set for IPAHFB=0. 

For IPABCS=1 the code HFODD solves the BCS equations with the neutron and proton pairing 
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strengths defined in Ref. [21J. These values can be modified by defining the multiplicative 
factors FACTGN and FACTGP for neutrons and protons, respectively. For IPABCS=2, the BCS 
pairing calculations are performed with fixed values of the neutron or proton pairing gaps 
equal to DELFIN or DELFIP, respectively For IPABCS=3, the BCS pairing calculations are per- 
formed with state-dependent pairing gaps corresponding to the pairing matrix elements (see 
the keyword PAIFLMATRI) calculated for the contact forces defined in the same way as for the 
HFB pairing calculations. 

Positive values of IPABCS require IPAIRI=1 and IPAHFB=0 and are incompatible with rotations 
IR0TAT=1 or broken simplex ISIMPY=0. IPABCS=2 requires IDEFIN=IDEFIP=1. 

Keyword: INI_INVERS 

0, = INIINV, INIKAR 

Allowed values of INIKAR=0, 1, 2, or 3 and INIINV=0, 1, 2, or 3 correspond to the following 
Dl transformations: 





INIKAR=0 INIKAR=1 INIKAR=2 INIKAR=3 


INIINV=0 
INIINV=1 
INIINV=2 
INIINV=3 


I R x R y R z 

P S X Sy S z 

T R x Ry R T Z 

pT cT cT cT 
r D x D y D z 



where P is the space inversion, T is the time reversal, Rk is the signature (rotation by angle 7r 
about the k = x, y, or z axis), and P T = PT, Sk = PRk (simplex), R£ = TRk (^-signature 7 ), 
and = TSk (A;-simplex T ). For INIKAR=INIINV=0, no transformation is performed and this 
option is inactive. 

Transformations are performed at the level of the densities, after the first iteration. As a secu- 
rity measure, nonzero values of INIKAR and INIINV require N0ITER=1. Such values also require 
SL0WEV=SL0W0D=SL0WPA=0.0, so as not to mix the old and new potentials corresponding to 
the original and transformed densities, respectively. They are also incompatible with IPRGCM^O. 

A given Dj h transformation must be accompanied by the correspondingly broken symmetries, 
that is, INIINV=1 or 3 requires IR0TAT=1 and INIINV=2 or 3 requires IPARTY=0. 

Keyword: INI_R0TAT 

0.0, 0.0, 0.0, = ALPINI, BETINI, GAMINI, INIR0T 
For INIR0T=1, the wave functions are rotated by the Euler angles a, (3, and 7 corresponding, 
respectively, to ALPINI, BETINI, and GAMINI (all in degrees). Transformations are performed 
at the level of the densities, after the first iteration. As a security measure, INIR0T=1 requires 
N0ITER=1. It also requires SL0WEV=SL0W0D=SL0WPA=0.0, so as not to mix the old and new 
potentials corresponding to the original and transformed densities, respectively. INIR0T=1 is 
incompatible with IPRGCM^O and requires a spherical HO basis of %u x = huj y = hu z . 

A rotation about the z axis must be accompanied by the broken signature, that is, ALPINI^O 
or GAMINI^O requires ISIGNY=0. 
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Keyword: PROJECTGCM 

0, 0, 0, 1, 1, 0, 1, 1, 

IPRGCM, IPROMI, IPROMA, NUAKNO, NUBKNO, KPROJE, 

IFRWAV, ITOWAV, IWRWAV 

For IPRGCM=1 and 2, the code calculates the diagonal and non-diagonal GCM kernels, respec- 
tively. In addition, for NUAKNO^ 1 or NUBKNO^ 1, the AMP kernels are calculated, as described 
in Section 12.11 although at present this option is not yet available for non-diagonal kernels 
of IPRGCM=2. The AMP is performed for doubled angular momenta from IPROMI to IPROMA 
and requires a spherical HO basis of huj x = tkv y = hu z . For even (odd) particle number 
IN_FIX+IZ_FIX, IPROMI, IPROMA, and KPROJE must be even (odd). For IPRGCM=1 and 2, the 
present version (v2.38j) also requires IR0TAT=1 and IPAIRI=0. 

The AMP has been tested up to the values of angular momenta of 70h, and therefore, IPROMA 
must not be larger than 2*70=140. After further tests, higher values could be used by increas- 
ing the parameter JMAX=70 in function dsmalg, which calculates the Wigner functions, see 
Eqs. (pQ), (H, and ©. 

NUAKNO is the number of G-T nodes, which are used to perform integrations over the a and 7 
Euler angles. For NUAKN0=1, these integrations are not performed (ID AMP) and the states are 
assumed to be axial with the doubled projection of the angular momentum on the z axis equal 
to KPROJE. NUBKNO is the number of G-L nodes, which are used to perform integrations over 
the P Euler angle. For NUAKN0>1 and NUBKN0>1, a full 3D AMP is performed and the value of 
KPROJE is ignored. NUAKN0>1 requires ISIMPY=0 and ISIGNY=0. IPROMA must be larger than 
the absolute value of KPROJE. 

For IPRGCM=2, the code calculates the GCM kernels between the states labeled by three-digit 
indices from "000" to "999". Indices of the "left" states vary between IFRWAV and ITOWAV, 
and these states are read from the disc. The index of the "right" state equals to ITOWAV, and 
this state is equal to the current state. In addition, for IWRWAV=1, the current state is saved 
on the disc with the index of ITOWAV. This allows for a simultaneous calculation of the "right" 
state along with calculating its kernels with all previously calculated and stored "left" states. 
The states can also be stored on disc without calculating kernels in the given run, that is, by 
setting the value of ITOWAV along with IWRWAV=1 and IPRGCM=0. Names of files on the disc 
are composed by concatenating the three-digit index, and FILWAV. 

Keyword: SAVEKERNEL 

= ISAKER 

For ISAKER=1, the code attempts reading the kernel files Nxxx-Lyyy-Rzzz-/ /FILKER, where the 
three-digit indices are: 

• xxx is the consecutive index of the kernel file, 

• yyy is the index of the left wave function, 

• zzz is the index of the right wave function. 
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In the work directory, the file names for all indices xxx are scanned, starting from 001. The 
kernels stored in these files are read into memory and are not recalculated. The kernels that 
have not been found in the kernel files are calculated and stored in the kernel file with the 
lowest available index xxx. In this way, one can submit many parallel jobs, see the keyword 
PARAKERNEL, that calculate kernels for different values of the Euler angles a, (3, and 7. The 
results are then collected in different kernel files with indices xxx attributed automatically. If 
any of the jobs is terminated before completing its task, the same input data can be resubmitted 
and the calculation automatically continues from the point where it has been interrupted. Once 
all the kernels will have been calculated, which requires a large CPU time, the AMP can be 
performed within a very small CPU time by reading, again automatically, all the created kernel 
files. ISAKER=1 requires IPRGCM>0 

Keyword: CHECKERNEL 

1 = ICHKER 

The names of kernel files are saved within these files. As a security measure, when reading the 
kernel files, their names are cross-checked against the saved information. This cross-checking 
can be deactivated by using ICHKER=0. This option is useful whenever the kernel files have 
been renamed for any reason. 

Keyword: PARAKERNEL 

0, 1, 1, 1, 1 = IPAKER, NUASTA, NUASTO, NUGSTA, NUGSTO 
For IPAKER=1, the code only calculates kernels for different values of the Euler angles a, (3, and 
7 and the AMP is suspended. Calculations are performed for the G-T nodes in the Euler angle 
a from NUASTA to NUASTO, for those in the Euler angle 7 from NUGSTA to NUGSTO, and for all the 
G-L nodes in the Euler angle f3, that is, from 1 to NUBKNO. For IPAKER=1, to save memory the 
code can be compiled with IPARAL=1. IPAKER=1 requires IPRGCM>0 and ISAKER=1. Values 
of NUASTA, NUASTO, NUGSTA, NUGSTO must all be between 1 and NUAKNO and must be ordered as 
NUASTA< NUASTO and NUGSTA<NUGST0. 

Keyword: TRANSITION 

2, 1, = NMURED, NMARED, NSIRED 
Maximum numbers of transition electric, magnetic, and surface or Schiff moments, respectively, 
for which proton kernels and reduced matrix elements are calculated, printed, and stored in 
the kernel files. NMARED and NSIRED must not exceed NMURED. For NMURED=0, NMARED=0, or 
NSIRED, the corresponding proton kernels and reduced matrix elements are not calculated. 

Keyword: CUTOVERLAP 

0, 10" 10 , 1. = ICUTOV, CUTOVE, CUTOVF 
For ICUT0V=0, parameters CUTOVE and CUTOVF are ignored and the collective states for the 
f^-mixing calculation, Eq. (fTO]) . are selected by their norm eigenvalues n m being larger than 
the negative of the smallest norm eigenvalue OVEMIN. For ICUT0V=1, the collective states are 
selected by their norm eigenvalues being larger than CUT0VE+CUT0VF*0VEMIN. 

Keyword: LIPKIN 

0, = LIPKIN, LIPKIP 

For LIPKIN=1 and/or LIPKIP=1, the Lipkin-Nogami corrections are included for neutrons 
and/or protons, respectively, see Section 12.91 At present, LIPKIN=1 or LIPKIP=1 requires 
IPAHFB=1. 
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Keyword: INI_LIPKIN 

0.1, 0.1 = FE2INK0), FE2INK1) 
For IC0NTI=0 or IC0NTI=1 and ILC0NT=0, the Lipkin-Nogami calculations are started by using 
the initial values of the neutron and proton Lipkin-Nogami parameters A2 fl98|) . FE2INI (0) and 
FE2INK1), respectively. For IC0NTI=0 or ILC0NT=0, values of FE2INI are ignored. 

Keyword: FIXLAMB2_N 

0.1, = FE2FIN, IF2FIN 
For IF2FIN=1, the neutron Lipkin-Nogami calculations are performed by using the fixed value of 
the neutron Lipkin-Nogami parameter A2 fl98|) equal to FE2FIN. IF2FIN=1 requires LIPKIN=1. 

Keyword: FIXLAMB2_P 

0.1, = FE2FIP, IF2FIP 
Same as above but for the proton Lipkin-Nogami calculations. 

Keyword: SL0WLIPKIN 

0. 5 = SL0WLI 

SL0WLI gives the value of the mixing fraction e used for the Lipkin-Nogami parameter A2 fl98|) . 
in analogy to the SL0WEV, SL0W0D and SL0WPA parameters. 

Keyword: FIXFERMI_N 

-8.0, = FERFIN, IFEFIN 
For IFEFIN=1, the HFB pairing calculations are performed with a fixed value of the neutron 
Fermi energy equal to FERFIN. At present, IFEFIN=1 requires IPAHFB=1. 

Keyword: FIXFERMI_P 

-8.0, = FERFIP, IFEFIP 

Same as above but for the proton HFB pairing calculations. 

3.3 Configurations 
Keyword: BL0CKSIZ_N 

1, = INSIZN, IDSIZN 

For |IDSIZN|=1, the code performs the neutron quasiparticle blocking calculations in the case 
when no symmetries are conserved, see Section [2771 For IDSIZN=+1 or —1, the blocked quasi- 
particle state is selected by having the largest overlap with the INSIZN-th neutron s.p. eigenstate 
of the HFB mean-field Routhian or with its time-reversed partner, respectively. Note that for 
rotating states, the time-reversed eigenstate is not necessarily an eigenstate of the Routhian. 
|IDSIZN| = 1 requires ISIMPY=0, IPARTY=0, IPAHFB=1, and IR0TAT=1. 

Keyword: BL0CKSIZ_P 

1, = INSIZP, IDSIZP 

Same as above but for the proton quasiparticle blocking. For odd-odd nuclei, neutron and 
proton quasiparticles can be simultaneously blocked. 

Keyword: BL0CKSI1VLN 

1, 0, = INSIMN, IRSIMN, IDSIMN 
For |IDSIMN|=1, the code performs the neutron quasiparticle blocking calculations in the case 
when simplex is conserved, see Section 12.71 For IDSIMN=+1 or —1, the blocked quasiparti- 
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cle state is selected by having the largest overlap with the INSIMN-th neutron s.p. eigenstate 
of the HFB mean-field Routhian in a given simplex block or with its time-reversed partner, 
respectively. The simplex of the state, +i or —i, is defined by IRSIMN=0 or 1, respectively. 
|IDSIMN| = 1 requires ISIMPY=1, IPARTY=0, IPAHFB=1, and IR0TAT=1. 

Keyword: BLOCKSIM_P 

1, 0, = INSIMP, IRSIMP, IDSIMP 

Same as above but for the proton quasiparticle blocking. For odd-odd nuclei, neutron and 
proton quasiparticles can be simultaneously blocked. 

Keyword: BL0CKSIQ_N 

1, 0, = INSIQN, IPSIQN, IDSIQN 
For |IDSIQN|=1, the code performs the neutron quasiparticle blocking calculations in the case 
when parity is conserved. For IDSIQN=+1 or —1, the blocked quasiparticle state is selected by 
having the largest overlap with the INSIQN-th neutron s.p. eigenstate of the HFB mean-field 
Routhian in a given parity block or with its time-reversed partner, respectively. The parity of 
the state, +1 or —1, is defined by IPSIQN=0 or 1, respectively. |IDSIQN| = 1 requires ISIMPY=0, 
IPARTY=1, IPAHFB=1, and IR0TAT=1. 

Keyword: BL0CKSIQ_P 

1, 0, = INSIQP, IPSIQP, IDSIQP 

Same as above but for the proton quasiparticle blocking. For odd-odd nuclei, neutron and 
proton quasiparticles can be simultaneously blocked. 

Keyword: BL0CKSIG_N 

1, 0, 0, = INSIGN, IPSIGN, ISSIGN, IDSIGN 
For |IDSIGN|=1, the code performs the neutron quasiparticle blocking calculations in the case 
when parity and signature are conserved. For IDSIGN=+1 or —1, the blocked quasiparticle 
state is selected by having the largest overlap with the INSIGN-th neutron s.p. eigenstate of the 
HFB mean-field Routhian in a given parity-signature block or with its time-reversed partner, 
respectively. The parity of the state, +1 or —1, is defined by IPSIGN=0 or 1, respectively. The 
signature of the state, +i or —i, is defined by ISSIGN=0 or 1, respectively. |IDSIGN| = 1 requires 
ISIMPY=1, IPARTY=1, IPAHFB=1, and IR0TAT=1. 

Keyword: BL0CKSIG_P 

1, 0, 0, = INSIGP, IPSIGP, ISSIGP, IDSIGP 

Same as above but for the proton quasiparticle blocking. For odd-odd nuclei, neutron and 
proton quasiparticles can be simultaneously blocked. 

Keyword: BL0CKFIX_N 

0, = IFIBLN, INIBLN 

For IFIBLN=1, the neutron quasiparticle blocking is based on calculating overlaps with a fixed 
s.p. wave function. The method is based on beginning the iteration by defining the number of 
the s.p. state, that is, INSIGN, INSIMN, INSIQN, or INSIZN, depending on the selected symmetry 
conditions. Then, for INIBLN=1 this one s.p. wave function is stored in memory and on the 
record file, and in consecutive iterations the overlaps are calculated with respect to this fixed s.p. 
wave function. Therefore, subsequent changes in the ordering and structure of s.p. states do not 
affect the blocking mechanism. For INIBLN=0, in the first iteration this fixed s.p. wave function 
is not initialized, but it is read from the record file. IFIBLN=1 requires |IDSIGN| = 1, |IDSIMN| = 1, 
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|IDSIQN|=1, or |IDSIZN|=1, depending on the selected symmetry conditions. IFIBLN=1 and 
INIBLN=0 requires IC0NTI=1. 

Keyword: BLOCKFIX_P 

0, = IFIBLP, INIBLP 

Same as above but for the proton quasiparticle blocking. 

3.4 Miscellaneous parameters 
Keyword: BROYDEN 

0, 7, 0.8, 1000. = IBR0YD, N_ITER, ALPHAM, BR0TRI 
For IBR0YD=1, the Broyden method is used to accelerate the convergence, see Section [2T8l For 
N_ITER=0, variables N_ITER, ALPHAM, and BR0TRI, which are read on this line, are ignored, that 
is, the previous values are kept. N_ITER is the number of iterations used to approximate the 
inverse Jacobian, ALPHAM is the value of the parameter a of the linear mixing to which is added 
the Broyden correction, and BR0TRI triggers an automatic switch to the Broyden method, that 
is, when the absolute value of the stability energy becomes lower than BR0TRI, iterations are 
changed to the Broyden scheme. A large value of BR0TRI ensures that the Broyden method 
is used from the very first iteration. IBR0YD=1 is incompatible with I_YUKA>0, IC0UDI=2, or 
IC0UEX=2. The Broyden method is implemented only in the FORTRAN-90 version of hfodd. 

3.5 Constraints 
Keyword: MULTC0NSCA 

1, 0., 0., = LAMBDA, STIFFG, GASKED, IFLAGG 

For IFLAGG=1, the total scalar multipole moment of the given multipolarity A is constrained. 
Values of LAMBDA, STIFFG, and GASKED correspond, respectively, to A, C\, and Q\ in Eq. 
(I80p . For IFLAGG=0, there is no scalar constraint in the given multipolarity. The constrained 
multipolarity LAMBDA cannot be equal to or exceed NMUC0N. For conserved parity IPARTY=1, 
only even moments can be constrained. 

Keyword: N0RBC0NSTR 

= N0_0RB 

For N0_0RB=1, constraints on the intrinsic spin only are used and the orbital part of the angular 
momentum is not constrained. N0_0RB=1 requires IR0TAT=1. 

Keyword: B0HR_BETAS 

4, 0, 1 = NEXBET, IPRIBE, IPRIBL 

For a given set of electric multipole moments, the code calculates and prints the corresponding 
first-order (for IPRIBL=1) and/or exact (for IPRIBE=1) Bohr deformation parameters, see 
Section I2~5l For IPRIBE=1, in case the code fails to find the exact values, the first-order values 
are printed irrespective of the value of IPRIBL. The exact values are sought for multipolarities 
up to NEXBET, which must not be greater than NMUPRI. The approximate values are printed up 
to multipolarity of NMUPRI. 

Keyword: SCHIFF_M0M 

= ISCHIF 
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For ISCHIF=1, the surface moments are everywhere in the code replaced by the Schiff moments, 
see Section [221 

3.6 Output-file parameters 
Keyword: HFBMEANFLD 

= IMFHFB 

For IMFHFB=1, eigenvalues of the HFB mean-field s.p. Hamiltonian or Routhian are printed. 
IMFHFB=1 requires IPAHFB=1. 

Keyword: PRINTJVMP 

0, 999, 1, 1, 1, 0, 

ISLPRI, ISUPRI, IENPRI, ISRPRI, IMIPRI, IKEPRI, IRMPRI 

These parameters govern the printing of the AMP results. If not restricted by the values of 
the doubled angular momenta for which the calculations are performed, IPROMI and IPROMA, 
the results are printed only for the doubled angular momenta between ISLPRI and ISUPRI. For 
IENPRI=1, the AMP energies are printed. In addition, for IENPRI=2, the AMP kernels are also 
printed. For ISRPRI=1, the sum rules are printed and compared with the HF average values. 
For IMIPRI=1, the energies of the .ff -mixed states are printed. For IKEPRI=1 and/or IRMPRI=1, 
the proton reduced kernels and/or reduced matrix elements, respectively, are printed. 

Keyword: TRANCUTPRI 

0., 0., 0., = QMUCUT, QMACUT, QSICUT 
Values of the electric, magnetic, and surface or Schiff proton kernels and reduced matrix ele- 
ments are printed only if their absolute values are larger than, respectively, QMUCUT, QMACUT, 
QSICUT. This may avoid printing long lists of very small or zero values. 

Keyword: 0NE_LINE 

1 = I1LINE 

For IlLINE^ 0, a one-line convergence report is printed at each iteration. For I1LINE=1, the 
code prints the values of deformation 7, total angular momentum, total angular frequency u>, 
and angle between vectors of angular momentum and frequency. For I1LINE=2, the code prints 
the values of neutron and proton pairing gaps and Lipkin-Nogami parameters A2. 

Keyword: PRINT_VIOL 

= IVIPRI 

For IVIPRI=1, the code prints integrals of several symmetry-violating terms. IVIPRI=1 re- 
quires IR0TAT=1. 

Keyword: NILSS0NLAB 

3 = NILXYZ 

This option used for NILXYZ=1, 2, or 3, allows for printing the Nilsson labels defined with 
respect to the x, y, or z axis, respectively. This feature is useful when the symmetry axis of a 
nucleus and its spin are aligned along the x or y axis and not along the z axis, for which the 
standard Nilsson labels are defined. This is particularly important for analyzing configurations 
of band heads of odd nuclei within the conserved y-signature-symmetry limit. 
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3.7 Files 



Keyword: WAVEF_FILE 

HFODD.WFN = FILWAV 
CHARACTER*68 file name of the wave function file. Must start at the 13-th column of the data 
line. The binary wave function files must exist if IPRGCM=2, and will be read, see the keyword 
PROJECTGCM. 

Keyword: KERNELFILE 

HFODD.KER = FILKER 
CHARACTER*68 file name of the kernel file. Must start at the 13-th column of the data line. The 
binary kernel files are written and read if ISAKER=1, see the keyword SAVEKERNEL. 

Keyword: YUKAWASAVE 

1 = IWRIYU 

For IWRIYU=1 and I_YUKA>2, the Yukawa file is saved on disc after each iteration is completed. 
The file contains the matrix elements of the Yukawa mean field. For IWRIYU=0 and I_YUKA>2, 
the file is saved only once, after all iterations are completed. For IWRIYU=— 1, the file is 
never saved. The Yukawa calculations can be restarted by using the Yukawa file. Therefore, 
IYC0NT=1 requires setting IWRIYU=0 or 1 in the run that is to be continued. Equivalently, the 
Yukawa calculations can be restarted by using the fields file, which requires setting IWRIFI=0 
or 1 in the run that is to be continued. 

Keyword: REPYUKFILE 

HFODD.YUP = FILYUP 

CHARACTER*68 file name of the Yukawa file. Must start at the 13-th column of the data line. 
The binary Yukawa file with the name defined in FILREP must exist if IYC0NT=1, and will be 
read. If the filenames FILREP and FILREC are identical, the Yukawa file will be subsequently 
overwritten as a new Yukawa file. 

Keyword: RECYUKFILE 

HFODD.YUC = FILYUC 

CHARACTER*68 file name of the Yukawa file. Must start at the 13-th column of the data line. 
If IWRIYU=1, binary Yukawa file is written after each HF iteration. It contains complete 
information that allows restarting the Yukawa calculations in another run of the code. To 
restart, one has to specify IYC0NT=1 and provide the name of the file by defining FILREP. 

Keyword: LIPKINSAVE 

1 = IWRILI 

For IWRILI=1 and LIPKIN=1 or LIPKIP=1, the Lipkin-Nogami file is saved on disc after each 
iteration is completed. The file contains the matrix elements of the particle density matrices. 
For IWRILI=0 and LIPKIN=1 or LIPKIP=1, the file is saved only once, after all iterations are 
completed. For IWRILI=— 1, the file is never saved. The Lipkin-Nogami calculations can be 
restarted by using the Lipkin-Nogami file. Therefore, ILC0NT=1 requires setting IWRILI=0 
or 1 in the run that is to be continued. Equivalently, the Lipkin-Nogami calculations can be 
restarted by using the fields file, which requires setting IWRIFI=0 or 1 in the run that is to be 
continued. 
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Keyword: REPLIPFILE 

HFODD.LIP = FILLIP 

CHARACTER*68 file name of the Lipkin-Nogami file. Must start at the 13-th column of the data 
line. The binary Lipkin-Nogami file with the name defined in FILLIP must exist if ILC0NT=1, 
and will be read. If the filenames FILLIP and FILLIC are identical, the Lipkin-Nogami file will 
be subsequently overwritten as a new Lipkin-Nogami file. 

Keyword: RECLIPFILE 

HFODD.LIC = FILLIC 

CHARACTER*68 file name of the Lipkin-Nogami file. Must start at the 13-th column of the data 
line. If IWRILI=1, binary Lipkin-Nogami file is written after each HF iteration. It contains 
complete information that allows restarting the Lipkin-Nogami calculations in another run of 
the code. To restart, one has to specify ILC0NT=1 and provide the name of the file by defining 
FILLIP. 

Keyword: FIELD_SAVE 

-1 = IWRIFI 

For IWRIFI=1, the field file, is saved on disc after each iteration is completed. The file contains 
the matrix elements of the mean field. For IWRIFI=0, the file is saved only once, after all 
iterations are completed. For IWRIFI=— 1, the file is never saved. To restart calculations of the 
exact Coulomb exchange energy, the field file must exist. IFC0NT=1 requires setting IWRIFI=0 
or 1 in the run that is to be continued. 

Keyword: FIELDJDLD 

= IWRIOL 

For IWRI0L=1, in the last iteration the mean fields are not updated, that is, the mixing fractions 
e are set equal to 1, irrespective of values of the SLOWEV, SLOWOD and SLOWPA parameters. In this 
way, information stored on the record file corresponds to the last but one iteration. Only then, 
restarting of the calculation from the field file leads to a smooth continuation of iterations. 

Keyword: REP_FIELDS 

HFODD.FIP = FILFIP 

CHARACTER*68 file name of the field file. Must start at the 13-th column of the data line. The 
binary field file with the name defined in FILFIP must exist if IFC0NT=1, and will be read. If 
the filenames FILFIP and FILFIC are identical, the field file will be subsequently overwritten 
as a new field file. 

Keyword: REC_FIELDS 

HFODD.FIC = FILFIC 

CHARACTER* 68 file name of the field file. Must start at the 13-th column of the data line. If 
IWRIFI=1, binary field file is written after each HF iteration. It contains complete information 
that allows restarting the calculations from the matrix elements of fields. To restart, one has 
to specify IFC0NT=1 and provide the name of the file by defining FILFIP. 
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3.8 Starting the iteration 
Keyword: CONTYUKAWA 

= IYCONT 

For IYC0NT=1, results stored in the Yukawa file are used to define the Yukawa fields in the 
first iteration; otherwise the Yukawa fields are set equal to zero. When the Yukawa fields are 
taken into account (i.e., for I_YUKA=2), and if a smooth restart and continuation of iterations 
from previously stored results is required, value of IYC0NT=1 must be used. IYC0NT=1 is 
incompatible with either of I_YUKA<2 or IC0NTI=0. 

Keyword: CONTLIPKIN 

= ILCONT 

For ILC0NT=1, results stored in the Lipkin-Nogami file are used to define the density matrix 
required for the Lipkin-Nogami calculations in the first iteration; otherwise the required den- 
sity matrix is set equal to zero. For the Lipkin-Nogami calculations (i.e., for LIPKIN=1 or 
LIPKIP=1), and if a smooth restart and continuation of iterations from previously stored re- 
sults is required, value of ILC0NT=1 must be used. ILC0NT=1 is incompatible with either of 
LIPKIN=LIPKIP=0 or IC0NTI=0. 

Keyword: CONTFIELDS 

= IFCONT 

For IFC0NT=1, results stored in the field file are used to define the matrix elements of fields in 
the first iteration; otherwise the matrix elements are recalculated. For the Gaussian-expansion 
method used to calculate the Coulomb energy and Coulomb mean field, that is for IC0UDI=2 or 
IC0UEX=2, and if a smooth restart and continuation of iterations from previously stored results 
is required, value of IFC0NT=1 must be used. IFC0NT=1 is incompatible with IC0NTI=0. 

4 Output file 

Together with the FORTRAN source code in the file hfodd.f , two examples of the output file 
are provided in ge064-a.out and snl20-b.out. Selected lines from the file ge064-a.out are 
presented in section TEST RUN OUTPUT below. This output file corresponds to the input 
file ge064-a.dat, reproduced in section TEST RUN INPUT below. 
The output file ge064-a.out contains the following new sections: 

Section CODE COMPILED WITH THE... lists the values of the most important array dimensions 
and switches declared in the PARAMETER statements, fixed at the compilation stage. 

Sections BOHR DEFORMATIONS (EXACT MULTIPOLE MOMENTS) give values of the Bohr deforma- 
tion parameters, determined as described in Section 12.51 

Section KERNELS AND AVERAGE VALUES... gives values of the diagonal norm kernels N KK 
and average AMP energies Hkk/Nkk for different values of I and K. Real and imaginary 
parts of N KK are both printed, although all imaginary parts should be equal to zero when the 
numerical precision is sufficiently good. 



36 



Section SUM RULES... compares real and imaginary parts of the sum rules (149]) and (I50j) with the 
corresponding HF average values. Results are printed for the norm kernels which must 

sum up to 1. This condition is a primary test of whether a sufficient number of the angular 
momenta are included in the sum of Eq. (jUJ). In the example of the file ge064-a. out, the sum 
rule for the norm kernel N0RM= . 3248 indicates a much too small value of I max and too small 
numbers of the G-T and G-L integration nodes. 

Sections REDUCED KERNELS... give values of the reduced kernels ( 1561) of electric and magnetic 
transition operators. 

Sections RESULTS OF THE K-MIXING... give values of the energies Ei of i^-mixed states, see 
Eq. and values of norm eigenvalues n m , see Eq. (jUJ). 

Sections REDUCED MAT.ELEMS... give values of the reduced matrix elements floTj) of electric and 
magnetic transition operators calculated for the i^-mixed states. 

Section NUMBERS OF CALLS TO SUBROUTINES gives the statistics of calls to subroutine, which 
together with the section EXECUTION TIMES IN SUBROUTINES illustrates the work flow of the 
code. 

5 FORTRAN source file 

The FORTRAN source code is provided in the file hfodd.f and can be modified in several 
places, as described in this section. 

5.1 FORTRAN-90 version 

Similarly as for the previous version (v2.08k), the code hfodd version (v2.38j) is written in 
FORTRAN- 77 and FORTRAN-90. In the FORTRAN source code provided in the file hfodd . f , 
all the FORTRAN-90 features are commented out and inactive. However, very simple modifica- 
tions of the source code can easily be performed to transform the code hfodd to FORTRAN-90, 
as described in Section (IV-5.2). 

The present version of the code hfodd version (v2.38j) is the last one that works under 
FORTRAN- 77; future releases will only use FORTRAN-90 programming. 

A set of c-shell and ex-editor scripts is provided within the hfodd distribution file, which 
allows an easy installation, compilation, and execution of the FORTRAN-90 version on a Linux 
computer running the INTEL© FORTRAN COMPILER. 

5.2 Library subroutines 

5.2.1 BLAS. 

The code hfodd requires an implementation of the BLAS (Basic Linear Algebra Sub- 
progams) interface for common dense vector and matrix computations. A reference implemen- 
tation of these subroutines can be downloaded from 
|http:/ /www.netlib.org/blas/blas.tgz| 

but for peak performance it is recommended that an optimized version of these routines should 
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be installed. Optimized machine-specific BLAS (such as Sun Performance Library and Intel 
Math Kernel Library) are available from most hardware vendors, and there are also third-party 
optimized implementations such as GOTO and ATLAS available for various architectures. 

The BLAS subroutines are in the REAL*8/C0MPLEX*16 version, and should be compiled 
without promoting real numbers to the double precision. On the other hand, the code hfodd 
itself does require compilation with an option promoting to double precision. Therefore, the 
code and the BLAS package should be compiled separately, and then should be linked together. 

5.2.2 Diagonalization subroutines. 

The code hfodd requires an external subroutine that diagonalizes complex hermitian ma- 
trices. Version (vl.60r) (see II) has been prepared with an interface to the NAGLIB subroutine 
f02axe, version (vl.75r) (see III) with an interface to the LAPACK subroutine zhpev, and 
version (v2.08i) (see IV) with an interface to the LAPACK subroutine zhpevx. In the present 
version (v2.38j), all these interfaces remain supported and can be activated as described in 
ILTV However, the recommended interface is now to the LAPACK subroutine zheevr, as 
described in this section. 

In version (v2.38j) we have implemented interface to the LAPACK subroutine zheevr, 
which can be downloaded (with dependencies) from 

http : //www. net lib . org/ cgi-bin/netlibf iles .pl?f ilename=/lapack/ complexl6/zheevr . f 
This subroutine works with unpacked matrices and hence performs calculations in less CPU 
time at the expense of using a larger memory. Alternatively, the entire LAPACK library can 
be downloaded from 

http : //www.netlib . org/lapack/lapack. tgz 

Subroutine ZHEEVR and its dependencies are in the REAL*8/C0MPLEX*16 version, and should 
be compiled without promoting real numbers to the double precision. Therefore, the code and 
the zheevr package should be compiled separately, and then should be linked together. 

In order to activate the interface to the LAPACK zheevr subroutine, the following modi- 
fications of the code hfodd (v2.38j) have to be made: 

1. Change everywhere the value of parameter I_CRAY=1 into I_CRAY=0. 

2. Change everywhere the value of parameter IZHPEV=0 into IZHPEV=3. 

3. If your compiler and linker do not support undefined externals, or subroutines called with 
different parameters, remove calls to subroutines cgemm, f02axe, zhpev, and zhpevx. 

5.2.3 Matrix inversion subroutines. 

The code hfodd requires external subroutines that invert complex matrices and calculate 
their determinants. In version (v2.38j) we have implemented interfaces to the LINPACK sub- 
routines zgedi and zgeco, which can be downloaded from 
http://www.netlib.Org/linpack/zgedi.f , and 

http: / / www.netnb.org/linpack/zgecolfl , respectively, together with 
http : / / www . net lib . org/linpack / zgef a. f | 

These subroutines are in the REAL*8/C0MPLEX*16 version, and should be compiled without 
promoting real numbers to the double precision. Therefore, the code and these subroutines 
should be compiled separately, and then should be linked together. 
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TEST RUN INPUT 



I This file (ge064-a.dat) contains input data for the code HFDDD (v2.38j) I 



NUCLIDE 

ITERATIONS 

ITERAT_EPS 

MAXANTIOSC 

PING-PONG 

CHAOTIC 

PHASESPACE 

SKYRME-SET 
SKYRME-STD 

SIMPLEXY 

SIGNATUREY 

PARITY 

ROTATION 

TSIMPLEX3D 

PAIRING 

HFB 

VACSIG_NEU 
VACSIG_PRO 
OPTI_GAUSS 

BASIS_SIZE 
SURFAC_PAR 
SURFAC_DEF 
SURFAC_DEF 

OMEGAY 

MULTCONSTR 

MULTCONSTR 



General data 



32 32 
20 

0.000001 
5 

0.0 3 



SIII 

1 1 

1 
1 

-1 



1 1 










Interaction 



Symmetries 



PPSP PPSM PMSP PMSM 

7 7 9 9 

PPSP PPSM PMSP PMSM 

9 9 



Configurations 





7 


7 


1 






14 


680 


800. 


32 


32 


1.23 


2 





0.00 


4 





0.00 



Constraints 



0.00 
2 



PRINT-ITER 










1 





1 


PRINT-MOME 
















1 


PRINT-INTR 









EALLMINMAX 


-40. 


0. 




EQUASI_MAX 


10.0 






MAX_MULTIP 










2 


4 


6 


BOHR_BETAS 










6 


1 


1 



1.0 2.7 1 

2 1.0 -1.3 1 

Output-file parameters 



REVIEW 
RECORDFILE 

RECORDSAVE 
REPLAYFILE 

COULOMFILE 



Files 





ge064-a.rec 
1 

ge064-a.rec 
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COULOMSAVE 
RESTART 

EXECUTE 

ITERATIONS 

MULTCONSTR 

MULTCONSTR 

BROYDEN 

RESTART 

EXECUTE 

ROTATION 

TSIMPLEX3D 

OMEGAY 

EXECUTE 

ITERATIONS 

SIMPLEXY 

SIGNATUREY 

TSIMPLEX3D 

PROJECTGCM 

SAVEKERNEL 

PRINT_AMP 

TRANCUTPRI 

CUTOVERLAP 

RECORDSAVE 

KERNELFILE 

BROYDEN 
EXECUTE 
ALL_DONE 



ge064-a. cou 
1 1 



5000 



0.575 




10 



Starting the iteration 





10 10 



10 1 1 
0.001 0.001 
1 0.01 
1 

ge064-a.ker 
0.0 



1 1 
0.001 
0.00 



0.0 



Calculate 
Next Run 



2 





1. 


.0 


2. 


.7 





2 


2 


1. 


,0 


-1. 


.3 





1 





0, 


.0 


0. 


.0 




1 














1 














1 


-1 




1 









Next Run 



Next Run 



Terminate 



42 



TEST RUN OUTPUT 



* HFDDD HFODD HFODD HFODD HFDDD HFDDD HFDDD HFDDD * 

* SKYRME-HARTREE-FOCK-BDGDLYUBOV CODE VERSION: 2.38J * 

* NO SYMMETRY-PLANES AND NO TIME-REVERSAL SYMMETRY * 

* DEFORMED CARTESIAN HARMONIC-OSCILLATOR BASIS * 

* * 

* * 

* J. DOBACZEWSKI, B.G. CARLSSON, J. DUDEK, J. ENGEL * 

* P. OLBRATOWSKI , P. POWALOWSKI , M. SADZIAK, W. SATULA * 

* N. SCHUNCK, A. STASZCZAK, M. STOITSOV, M. ZALEWSKI * 

* AND H. ZDUNCZUK * 

* * 

* INSTYTUT FIZYKI TEORETYCZNEJ , WARSZAWA * 

* 1993-2009 * 
****************** ********** 

* * 

* CODE COMPILED WITH THE FOLLOWING ARRAY DIMENSIONS AND SWITCHES: * 



NDBASE = 

NDMAIN = 
NDAKNO = 
IPARAL = 



680 
14 
50 




NDSTAT 
NDMULT 
NDBKNO 
I CRAY 



181 
9 
50 




NDXHRM = 
NDMULR = 
NDPROI = 
IZHPEV = 



30 
4 

50 
3 



NDYHRM = 30 NDZHRM = 30 
NDLAMB = 9 NDITER = 5000 
NDCOUL = 80 NDPOLS = 25 



******************************************************************************* 

* REDUCED KERNELS AND MATRIX ELEMENTS: * 

* FOR ELECTRIC MULTIPOLES CALCULATED UP TO = 2, PRINTED IF LARGER THAN 0.001 * 

* FOR MAGNETIC MULTIPOLES CALCULATED UP TO = 1, PRINTED IF LARGER THAN 0.001 * 

* FOR SURFACE MULTIPOLES CALCULATED UP TO = 0, PRINTED IF LARGER THAN 0.001 * 

*** ****** 



* PRINTING THE RESULTS FOR ANGULAR-MOMENTUM-PROJECTED STATES: < 2*1 < 10 



AVERAGE ENERGIES OF PROJECTED STATES 
KERNELS BETWEEN THE PROJECTED STATES 
SUM RULES COMPARED TO THE H-F VALUES 
K-MIXED ENERGIES OF PROJECTED STATES 
REDUCED KERNELS BETWEEN THE PROJECTED STATES 
REDUCED MATRIX ELEMENTS BETWEEN THE K-MIXED STATES 



YES 
NO 
YES 
YES 
YES 
YES 



* BOHR DEFORMATIONS (EXACT MULTIPOLE MOMENTS) TOTAL * 

* * 



BIO = 
B20 = 
B30 = 
B40 = 
B50 = 

B60 = 



ZERO 
0.2094 
ZERO 
-0.0306 
ZERO 



Bll = 

B21 = 

B31 = 

B41 = 

B51 = 



ZERO 
ZERO 
ZERO 
ZERO 
ZERO 



B22 = 
B32 = 
B42 = 
B52 = 



-0.1138 
ZERO 

-0.0348 
ZERO 



-0.0054 B61 = 



ZERO B62 = 0.0147 



B33 = 
B43 = 
B53 = 

B63 = 
B65 = 



ZERO 
ZERO 
ZERO 

ZERO 
ZERO 



B44 = 
B54 = 
B55 = 
B64 = 



-0.0049 
ZERO 
ZERO 

-0.0098 



B66 =-7.7E-04 



* KERNELS AND AVERAGE VALUES OBTAINED FOR ANGULAR-MOMENTUM PROJECTED STATES * 

* * 

* * 

* REAL (KERNEL) 



IMAG (KERNEL) AVERAGE ENERGY 



* 


I,K= 








NORM = 





008035212882 





000000000000 


-542 


137728 


* 


I,K= 


1 


-1 


NORM = 





000867254021 





000000000000 


-509 


681812 


* 


I,K= 


1 





NORM = 





000000000000 





000000000000 


-967 


529412 


* 


I,K= 


1 


1 


NORM = 





000867254021 





000000000000 


-509 


681812 


* 


I,K= 


2 


-2 


NORM = 





014229645333 





000000000000 


-541 


456596 


* 


I,K= 


2 


-1 


NORM = 





003145537421 





000000000000 


-530 


120652 
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* I,K= 2 NORM = 0.034483057875 0.000000000000 -541.966975 * 

* I,K= 2 1 NORM = 0.003145537421 0.000000000000 -530.120652 * 

* I,K= 2 2 NORM = 0.014229645333 0.000000000000 -541.456596 * 

* SUM RULES OBTAINED FOR ANGULAR-MOMENTUM PROJECTED STATES VS THE HF VALUES * 

* * 



REAL (SUM RULE) IMAG(SUM RULE) 



NORM= 
SKYRME= 
EKIN_T= 
COUL_D= 
COUL_E= 

MULT IP OLE L= M= QMUL_P= 
L= 1 M=-l QMUL_P= 
L= 1 M= QMUL_P= 
L= 1 M= 1 QMUL_P= 
L= 2 M=-2 QMUL_P= 
L= 2 M=-l QMUL_P= 
L= 2 M= QMUL_P= 
L= 2 M= 1 QMUL_P= 
L= 2 M= 2 QMUL_P= 



0.324825071083 
-587.075736 
358.576775 
57.290170 
-4.008437 

10.394402 
0.000000 
0.000000 
0.000000 

-0.206203 
0.000000 
0.432164 
0.000000 

-0.206203 



0.000000000000 
0.000000 
0.000000 
0.000000 
0.000000 

0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.032421 
0.000000 
-0.032421 
0.000000 



HF VALUE 

-1808.113218 
1105.378182 
176.353367 
-12.345178 

32.000000 
0.000000 
0.000000 
0.000000 
0.000000 
0.000000 
1.333826 
0.000000 

-0.621629 



* REDUCED KERNELS OF PROTON MULT IP OLE MOMENTS [UNITS: (10 FERMI) "LAMBDA] * 

* ONLY THE RESULTS WITH ABSOLUTE VALUES LARGER THAN 0.001 ARE BELOW PRINTED * 

* * 

* * * 
* 

* RESULTS OF THE K-MIXING CALCULATION FOR THE FOLLOWING OVERLAP CUTOFF DATA * 

* ICUT0V=1, CUTOVE= 1.0E-02, CUTOVF= . OE+00 * 

* * 



IL 


KL 


L 


IR 


KR 


<l IQI l> 




IL 


KL 


L 


IR 


KR 


<l IQI l> 










2 


2 


-2 


-0.008 


* 


2 


-2 


2 


2 


-2 


0.025 


* 


2 





2 


2 


-2 


-0.017 


* 


2 


2 


2 


2 


-2 


0.025 


* 


1 


-1 


2 


2 


-1 


0.003 


* 


1 


1 


2 


2 


-1 


0.003 


* 








2 


2 





0.024 


* 


2 


-2 


2 


2 





-0.017 




2 





2 


2 





-0.055 




2 


2 


2 


2 





-0.017 


* 


1 


-1 


2 


2 


1 


-0.003 


* 


1 


1 


2 


2 


1 


-0.003 


* 








2 


2 


2 


-0.008 


* 


2 


-2 


2 


2 


2 


0.025 




2 





2 


2 


2 


-0.017 


* 


2 


2 


2 


2 


2 


0.025 


* 



I 


N 


OVERLAP 




ENERGY i 


I 


N 


OVERLAP 




ENERGY 


2 


1 


4.02E-02 


-541 


817486 I 


3 


1 


3.78E-02 


-540 


997225 


2 


2 


2.71E-02 


-541 


344808 I 












4 


1 


7.90E-02 


-541 


775088 I 


5 


1 


5.06E-02 


-540 


367910 


4 


2 


3.04E-02 


-540 


196719 I 


5 


2 


2.04E-02 


-537 


665199 


4 


3 


2.06E-02 


-538 


486590 I 

i . 













* 
* 

* REDUCED MAT . ELEMS OF PROTON MULT IP OLE 

* ONLY THE RESULTS WITH ABSOLUTE VALUES 

* 

* 
* 

* NUMBERS OF CALLS TO SUBROUTINES 



-* 

MOMENTS [UNITS: (10 FERMI) "LAMBDA] * 
LARGER THAN 0.001 ARE BELOW PRINTED * 

* 
* 

* 
* 
* 



IL 

2 
2 
3 
2 
4 
2 
4 
4 
2 
4 
3 
5 
5 



NL 

1 
1 
1 
2 
2 
1 
1 
3 
2 
2 
1 
2 
1 



L IR 



2 
2 
3 
4 
4 
4 
4 
4 
4 
4 
5 
5 
5 



NR 


<l IQI l> 


* 


IL 


NL 


L 


IR 


NR 


<l IQI l> 


1 


-1.142 


* 


2 


2 


2 


2 


1 


-1.419 


2 


-1.419 


* 


2 


2 


2 


2 


2 


1.142 


1 


-0.007 


* 


2 


1 


2 


4 


1 


-2.353 


1 


0.543 


* 


4 


1 


2 


4 


1 


-0.518 


1 


-1.023 




4 


3 


2 


4 


1 


-0.417 


2 


-0.554 




2 


2 


2 


4 


2 


-1.264 


2 


-1.023 




4 


2 


2 


4 


2 


-2.302 


2 


0.763 


* 


2 


1 


2 


4 


3 


0.159 


3 


-0.831 


* 


4 


1 


2 


4 


3 


-0.417 


3 


0.763 




4 


3 


2 


4 


3 


2.832 


1 


2.170 


* 


5 


1 


2 


5 


1 


-0.937 


1 


1.296 


* 


3 


1 


2 


5 


2 


0.381 


2 


1.296 




5 


2 


2 


5 


2 


0.982 
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X** 

* 
* 

* 
* 

* 
* 

* 
* 



************* **** *************** Jit********************************** 



4001 => INTORB 2354 => BEGINT * 

2162 => DENMAC 2000 => RDTWAV * 

1162 => MAGMOM 1084 => CDUMAT * 

546 => INTSOR 541 => DIAMAT * 

249 => ROTMOM 162 => AVPARI * 

162 => LINAVR 160 => DIASIG * 

83 => RECORD 81 => SKFILD • 

81 => SHISIF 81 => SHIMAG * 

14 => AVOBSE 14 => AVSIMP * 

2 => NILSON 2 => DIASIZ * 

1 => PROANG 1 => DIAPRO * 

* 
* 
* 



26730 
2324 
1162 
1084 
358 
162 
160 
81 
59 
14 
1 
1 



=> INTMUL 
=> INTKIN 
=> MOMETS 
=> BEGINC 
=> INTMAS 
=> MOMVMU 
=> INTEGH 
=> SHICOE 
=> DOBROY 
=> NILABS 
=> HFODD 
=> REDPRO 



4002 
2162 
1162 
1084 
358 
162 
106 
81 
14 
5 
1 



INTSPI 
DENSHF 
MOMSIF 
INTCOU 
INTCEN 
SPIMOM 
SPAVER 
SHIMOM 



=> AVANGY 
=> RECOUL 
=> POWALL 
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